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INTRODUCTION  -  subject.  Traumatic  brain  injury  (TBI)  is  a  major  contributor  to  morbidity  and 
mortality  experienced  by  those  soldiers  subjected  to  improvised  explosive  devices  (IED)  as  well  as  in  military 
and  civilian  high  speed  collisions.  Traumatic  cerebral  vasospasm  (TCV)  is  a  major  contributor  to  morbidity  and 
mortality  experienced  by  those  TBI  patients.  Aggressive  neurosurgical  treatment  motivated  by  early  diagnosis 
appears  to  improve  the  clinical  outcome  for  these  patients.  Early  transcranial  Doppler  (TCD)  measurements  of 
blood  flow  speed  within  major  cerebral  arteries  produces  the  initial  diagnosis,  hence  motivates  the  rapid 
treatment  of  TCV.  However,  the  skill  necessary  to  deploy  TCD  limits  its  availability  relative  to  its  need.  This 
gap  in  patient  care  reduces  the  quality  of  care  and  potentially,  therefore,  the  quality  of  life  of  injured  soldiers. 
This  gap  also  defines  a  critical  need  for  robust,  easy  to  use  TCD  system,  applicable  to  these  patients. 

INTRODUCTION  -  purpose.  Early  work  by  the  PI  on  automating  transcranial  Doppler  has  been 
translated  by  PhysioSonics,  Inc,  (PSI)  a  Seattle-based  company  that  he  co-founded  into  a  working  automated 
transcranial  Doppler  (aTCD)  device  called  ‘Presto’  that  they  have  submitted  to  the  FDA  via  a  51  OK  for  approval. 
Their  system  is  based  upon  a  novel  and  proprietary  ultrasound  platform  along  with  novel  and  proprietary 
‘search  and  lock  on’  algorithms  for  detecting  and  staying  focused  on  the  brightest  spot  in  ‘Doppler’  space  -  that 
is,  within  ultrasound-derived  maps  of  blood  flow  speed  captured  by  their  device.  This  is  reasonable  because  a 
view  of  each  of  the  major  cerebral  arteries  relevant  to  vasospasm,  listed  below,  shows  that  each  contribute  the 
largest  values  of  blood  flow,  hence  the  largest  signal  in  Doppler  space.  PhysioSonics  has  optimized  it  to 
measure  blood  flow  speed  within  the  middle  cerebral  artery  (MCA)  -  sufficient  for  some  assays  of  vasospasm, 
and  demonstrated  its  utility  on  uninjured  volunteers.  However,  a  critical  additional  assay  of  vasospasm  -  the 
Lindegaard  ratio  -  requires  measurements  of  blood  flow  speed  within  the  extra-cranial  segment  of  the  internal 
carotid  artery  as  well  as  from  the  MCA.  Moreover,  while  their  system  already  works  on  humans,  PhysioSonics 
has  no  plans  to  optimize  the  device  for  TCV,  or  to  extract  information  from  the  ICA. 

INTRODUCTION  -  scope  of  the  research.  To  met  the  goal  of  this  proposal  we  sought  to  test  the 
following  hypothesis:  PhysioSonics’  automatic  TCD  (aTCD)  system,  with  modifications  to  its  hardware  and/or 
software  as  necessary,  will  match  measurements  of  cerebral  blood  flow  in  the  middle  cerebral  artery  (MCA) 
and  in  the  internal  carotid  artery  (ICA)  made  by  standard  TCD  on  patients  who  suffer  from  cerebral  vasospasm. 

Specific  Aim  #  1 :  To  test  the  ability  of  the  aTCD  device  to  measure  blood  flow  velocity  in  the  MCA  and 
ICA  of  patients  with  known  vasospasm,  as  compared  to  standard  TCD  run  by  a  trained  neurosonographer. 
Proposed  study:  Use  the  device  to  measure  blood  flow  in  the  MCA  and  ICA  of  civilian  vasospasm  patients  - 
typically  those  with  sub-arachnoid  hemorrhage  -  within  Harborview  Medical  Center  (HMC),  the  only  Level  One 
trauma  center  in  Seattle.  Compare  those  measurements  with  those  determined  during  routine  clinical  care  by 
the  resident  neurosonographer  with  their  standard  clinical  TCD  unit.  Determine  the  level  of  fidelity  of  the 
aTCD-derived  information  on  blood  flow  with  that  given  by  the  gold  standard,  operator-dependent  TCD  device. 

Specific  Aim  #2:  To  modify  as  necessary  the  TCD  array  and  associated  holder  for  use  on  the  ICA,  and, 
the  algorithms  required  to  extract  the  highest  blood  flow  velocities  in  the  MCA  and  ICA  to  guarantee  that  the 
aTCD  device  can  detect  the  highest  levels  of  vasospasm.  Proposed  study,  ill  Modify  as  necessary  the  brace 
that  ordinarily  holds  the  transducer  to  the  temple  (to  facilitate  measurement  of  blood-flow  from  within  the  MCA 
so  that  it  will  automatically  extract  blood  flow  information  from  the  ICA.  Proposed  study.  Modify  as 
necessary  the  software  within  the  aTCD  device  to  optimize  its  ability  to  calculate  the  value  of  fastest  blood  flow 
within  the  MCA  and  the  ICA  for  vasospasm  patients.  Proposed  study.  #3:  Modify  as  necessary  the  TCD 
transducer  or  its  controlling  software  so  as  to  optimize  its  ability  to  locate  the  point  of  fastest  blood  flow  within 
the  MCA  and  the  ICA  for  vasospasm  patients 

Impact.  Successful  completion  of  the  proposed  work  will  yield  a  device  that  could  readily  become 
available  to  the  military  with  its  data  interpreted  on  site  or  sent  or  via  telemedicine  for  interpretation  by  distant 
neurosurgical  experts. 
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BODY 


Specific  Aim  #1 :  To  test  the  ability  of  the  aTCD  device  to  measure  blood  flow  velocity  in  the  middle 
cerebral  artery  (MCA)  and  internal  carotid  artery  (ICA)  of  patients  with  known  vasospasm,  as  compared  to 
standard  TCD  run  by  a  trained  neurosonographer.  (Q1-Q10  out  of  12  quarters) 

Task  1 :  Get  approval  for  the  human  subjects  study  (Q1-2).  This  study  will  involve  up  to  128  patients 
with  vasospasm,  with  multiple  studies  per  patient,  over  the  three  years  of  the  study. 

(1.1)  Done  in  the  first  year. 

Task  2:  Take  possession  of  the  aTCD  system  (Q3). 

(2.1)  Done  in  the  first  year. 

Task  3:  Modify  the  software  associated  the  aTCD  system  so  it  will  report  and  store  the  spectrogram 
derived  from  each  of  the  MCA  and  ICA,  including  the  Lindegaard  ratio.  (Q2-5). 

(3. 1)  Done  in  the  first  year. 

Milestone  #1:  Achieve  ability  to  start  human  study.  Done. 

Task  4:  Modify  as  necessary  the  commercial  headband  that  holds  the  TCD  so  that  it  can  deploy  on  the 
upper  portion  of  the  neck  of  volunteers.  (Q1-2) 

(4. 1)  Done  in  the  first  year. 

Task  5:  Collect  then  archive  aTCD  and  standard  TCD  data  from  the  MCA  and  ICA  of  128  patients, 
along  with  appropriate  supporting  medical  data.  For  a  subset  of  these  patients  we  will  also  collect  standard 
TCD  data  twice  in  a  day  in  order  to  determine  the  length  of  time  we  will  need  to  eventually  monitor  patients  with 
our  aTCD  (Q3-Q10)  UW. 

(5. 1)  We  started  to  recruit  patients  with  civilian  vasospasm  in  the  5/1 2th  quarter  of  our  study  (October 
2012).  Table  1  documents  the  number  of  patients  identified,  consented,  and  studied  during  this  research  effort 
(October  201 2  -  August  201 5). 


Human  Data 

Identified 

Consented 

Studied 

#  of  data  sets 

1 0ct201 2- 
15Sept2013 

65 

22 

16 

25 

1Oct2013  - 
30Sept2014 

148 

37 

34 

46 

1Oct2014  - 
31  Aug2015 

125 

42 

33 

48 

Total 

338 

101 

83 

119 

Note  that  we  lose  two  thirds  of  our  identified  patients  to:  lack  of  bone  windows,  lack  of  available  next  of  kin, 
refusal  to  participate  (a  big  problem  because  these  patients  are  typically  awake  but  uncomfortable),  lack  of 
cranial  bone  and  patient  death. 
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Task  6:  Compare  spectrogram  features  found  within  aTCD  and  standard  TCD  data.  (Q5-Q10)  UW 
(6.1)  In  Figures  1-3  we  compare  representative  blood-flow  data  in  the  MCA  of  three  vasospasm 
patients.  (Note  that  we  get  pdfs  from  the  commercial  system  whose  resolution  for  this  report  fall  short  of 
adequate  representation  when  reproduced  here.)  These  examples  demonstrate  that  we  can  collect  TCD  data 
from  the  MCA  of  patients,  with  a  mix  of  success  in  their  comparison  with  the  gold  standard  data,  in  ways  that 
have  or  will  produce  improvements  in  the  aTCD  system. 

Our  system’s  data  falls  short  in  two  key  areas  -  the  search  algorithm  and  the  spectral  envelope 
calculation  -  as  anticipated  during  the  construction  of  this  grant  application.  Indeed,  our  tasks  as  originally 
designed  below  seek  to  address  these  issues.  As  described  in  our  2nd  annual  report,  we  modified  the  aTCD 
device  to  search  deeper,  where  our  previous  maximum  depth  was  50  mm  and  is  now  70mm.  We  have  now 
also  added  a  feature  where  the  user  is  able  to  force  the  machine  to  find  peak  flow  in  a  certain  volume  of  depth, 
by  moving  search  box  covering  10  mm  vertical  distance  -  this  is  described  in  detail  in  Task  10  below.  This  is 
meant  to  address  some  issues  of  the  locking  flow  selection  feature,  where  non-peak  flow  was  selected. 

Also,  here  we  show  the  standard  TCD  system’s  reports  of  mean  flow  speeds,  e.g.  107  cm/s  at  a  depth 
in  the  left  MCA  of  50  mm,  and  pulsatility  index  (PI),  a  relative  measure  of  the  tendency  of  the  peak  systolic 
blood  flow  speed  to  rise  above  the  mean  blood  flow  speed.  The  aTCD  system  also  reports  these  values,  but  in 
a  different  window.  We  will  develop  quantitative  comparisons  of  these  quantities  as  well  as  of  peak  systolic 
blood  flow  speed.  Also,  our  ICA  measurements  with  the  aTCD  system  compare  favorably  with  that  of  the 
clinical  system.  We  have  focused  on  those  measurements  that  motivate  changes  to  the  aTCD  system  here  for 
this  report. 


Figure  1.  Data  from  the  aTCD  machine  (A)  and  a 
standard  TCD  machine  (B),  where  the  aTCD  data 
under-represents  blood  flow  speeds  in  the  MCA. 

Here  in  (A)  we  see  peak  systolic  velocities  identified 
by  the  commercial  spectral  envelope  varying  between 
1 00-1 50  cm/s.  However,  in  (B)  the  figures  show  that 
the  standard  TCD  machine,  with  user  intervention, 
measures  peak  systolic  velocities  as  a  function  of 
distance  from  the  transducer  [at  50,  56  and  42  mm 
distance]  with  values  closer  to  180  cm/s.  The  aTCD 
machine  provides  data  suggesting  this  patient  is  in 
mild  vasospasm  while  the  clinical  machine 
demonstrates  that  this  patient  is  in  severe 
vasospasm.  Therefore,  for  this  patient,  the  flow- 
envelope  calculator  fails  to  capture  peak  flow  values. 
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Figure  2.  Data  from  the  aTCD  machine  (A)  and  a 
standard  TCD  machine  (B),  where  the  aTCD  data 
well  represents  the  blood  flow  speeds  in  the  MCA 
as  defined  by  the  standard  TCD  machine.  Here  in 
(A)  we  see  peak  systolic  velocities  identified  by  the 
commercial  spectral  envelope  (in  fucia)  140-150  cm/s, 
nominally  close  to  the  underlying  spectrogram  (in 
green).  This  compares  favorably  to  measurements  of 
peak  systolic  blood  flow  speeds  (B)  by  the  standard 
TCD  machine  as  a  function  of  depth  [at  40,  44,  50  and 
54  mm  distance  from  the  transducer].  These  peak 
speeds  start  at  a  value  of  120  cm/s  at  40  mm  distance 
from  the  transducer,  increasing  to  close  to  150  cm/s 
by  54  mm  distance.  Both  the  aTCD  machine  and 
clinical  machine  demonstrates  that  this  patient  is  in 
moderate  to  severe  vasospasm.  Therefore,  for  this 
patient,  the  aTCD  system  performs  plausibly  well 
compared  to  the  clinical  system. 
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Figure  3.  Data  from  the  aTCD  machine  (A)  and  a 
standard  TCD  machine  (B),  where  the  aTCD  data 
collects  representative  blood  flow  speeds  in  the 
MCA,  but  with  a  poor  spectral  envelope 
calculation,  and  at  too  shallow  a  depth,  as 
compared  to  the  data  collected  by  the  standard 
TCD  machine.  Here  in  (A)  we  see  peak  systolic 
velocities  identified  by  the  commercial  spectral 
envelope  (in  fucia)  110  cm/s,  woefully  low  compared 
to  the  150  cm/s  readily  identify  in  the  underlying 
spectrogram  (in  green).  These  ‘green’  values  of  blood 
flow  speed  compare  favorably  only  to  measurements 
of  peak  systolic  blood  flow  speeds  (B)  by  the  standard 
TCD  machine  at  its  shallower  depths  (between  44  and 
50  mm  distance  from  the  transducer  face).  The  aTCD 
search  algorithm  completely  misses  the  very  high  flow 
speeds  captured  by  the  standard  TCD  machine, 
namely,  at  or  over  240  cm/s  between  depths  of  50-64 
mm.  Therefore,  for  this  patient,  the  flow-envelope 
calculator  fails  to  capture  peak  flow  values  and  the 
device’s  search  algorithm  locks  onto  local  peak  flow 
velocities  that  are  too  shallow,  compared  to  the 
manually  steered  standard  TCD  system. 
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We  have  developed  three  enveloping  algorithms,  of  which  one  appears  to  track  best  with  standard 
clinical  TCD  machines.  Figure  4(a),  below,  shows  the  average  TCD  velocity  across  all  patients  compared 
between  the  Presto  algorithm,  our  G90  and  G95  algorithms,  and  compared  to  the  standard  clinical  machine 
data  (looking  at  shallow  flow  data  <  50  mm  deep,  and  at  the  maximum  across  the  entire  depth  range).  Figures 
4(B)  and  4(C)  are  two  examples  of  envelopes  calculated  using  our  three  new  algorithms  and  the  original  Presto 
algorithm.  The  yellow  envelope  (GMM  95/G95)  appears  to  track  best  with  standard  clinical  TCD  data,  when  we 
look  only  at  shallow  depths  (that  is,  less  than  50mm).  The  signal-locating  algorithm  on  the  Presto  has  a 
tendency  to  focus  on  signals  within  50mm  of  the  transducer  face,  therefore  we  believe  that  a  majority  of  our 
Presto  data  comes  from  within  that  range. 
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Figure  4(A).  Average 
maximum  velocity  across  all 
TCD  data. 

Presto  =  Presto  algorithm 
G90  &  G95  =  two  envelope 
algorithms  that  we  have 
developed 

TCD  Shallow  =  maximum 
velocity  for  standard  clinical 
TCD  at  depths  <50mm 
TCD  Deep  =  maximum 
velocity  for  standard  clinical 
TCD  at  depths  >50mm 
TCD  Max  =  maximum  velocity 
for  standard  clinical  TCD 
across  all  depths 


300 

250 

5T  200 


&  150 

o 

o 

<D 

>  100 


50 


0 


Presto  (53)  G90  (53)  G95  (53)  TCD  Shallow  (63)TCD  Deep  (63)  TCD  Max  (63) 


Presto 
GMM  90 
GMM  95 
GMM  99 


Figure  4(B).  Three  of  our  new 
envelope  algorithms  applied  to 
a  sample  data  set.  The  yellow 
envelope  (GMM  95)  appears 
to  track  best  with  clincal  TCD 
data.  Note  that  GMM95 
corresponds  to  G95  in  Fig 
4(A). 
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Figure  4(C).  Three  of  our  new 
envelope  algorithms  applied  to 
a  sample  data  set.  The  yellow 
envelope  (GMM  95)  appears 
to  track  best  with  clincal  TCD 
data.  Note  that  GMM95 
corresponds  to  G95  in  Fig 
4(A). 
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Figure  4(D).  An  example 
where  the  Presto’s  native 
envelope  algorithm  fails  to 
capture  the  spectrogram 
accurately.  The  G95  model 
(yellow  line)  is  much  more 
effective  and  accurate  in 
capturing  the  peaks  in  the 
data.  Note  that  GMM95 
corresponds  to  G95  in  Fig 
4(A). 
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We  also  analyzed  the  results  for  the  best  envelope  algorithm  compared  to  standard  TCD  results  as  a 
function  of  percent  difference  over  time  during  the  study.  The  results  show  a  steady  trend  toward  greater 
accuracy  in  our  Presto  results  compared  to  standard  TCD.  Our  initial  results  in  the  first  6  months  of  data 
collection  averaged  6.5%  off  from  the  standard  TCD.  By  18-24  months  into  the  study,  our  results  averaged  less 
than  1 .5%  off  from  the  standard  TCD  (Figure  5,  below). 
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Percent  Difference  Between  Absolute  MaxTCD  and  G95 
Values  (Number  of  Patients) 


Figure  5.  Percent  difference  of  our  Presto  results  vs  standard  clinical  TCD  data.  The  difference  decreases  over 
time  through  the  study  as  we  improved  the  Presto  device  and  our  spectral  envelope  calculatiions. 


Milestone  #2\  Develop  quantitative  comparisons  of  the  aTCD-derived  and  standard  TCD-derived  blood 
flow  data  from  vasospasm  patients,  thereby  highlighting  short  comings,  if  any,  of  the  aTCD  system  for  the  case 
of  vasospasm. 

We  have  collected  sufficient  number  of  data  sets  from  our  patients  to  identify  three  problems  with  our 
existing  aTCD  system:  [la]  the  flow-envelope  calculator  fails  to  capture  peak  flow  values  and  [1b]  the  device’s 
search  algorithm  locks  onto  local  peak  flow  velocities  at  too  shallow  a  depth  for  many  patients  with  vasospasm, 
all  relative  to  the  clinical  gold  standard.  In  addition,  thanks  to  the  formation  of  peripheral  edema  in  several  of 
these  patients,  we  can  say  that  [1c]  the  device’s  search  volume  is  too  shallow  for  an  appreciable  percentage  of 
patients  with  vasospasm.  Task  9  addresses  problem  #1  a  while  Task  10  addresses  problems  #1b,c. 

Task  7:  Write  up  results  for  publication  and  presentations.  (Q7-8;  Q10-12)  UW 

(7. 1)  The  Journal  of  Ultrasound  in  Medicine  has  published  our  paper  (Marzban  et  al,  2013)  that 
describes  our  first  robust  means  of  calculating  optimal  spectrogram  envelopes  for  different  patients.  We  have 
submitted  two  additional  papers.  One  describes  a  second  means  of  calculating  optimal  spectral  envelopes 
(Morison  et  al,  2015),  one  based  upon  a  first  principles  method  described  in  our  previous  report  and  the  results 
are  summarized  below  beginning  at  Figure  10.  This  paper  has  been  submitted  to  Biomedical  Signal 
Processing  and  Control.  We  have  published  a  third  paper  associated  with  this  project,  one  based  upon  a  new 
statistics-based  modeling  effort  for  calculating  spectra,  which  we  report  below  (Marzban  et  al,  2016). 

Milestone  #3:  get  feedback  on  the  quantitative  aspects  of  our  work  from  our  scientific  and  clinical  peers. 

Our  clinical  peers  are  quite  excited  by  the  prospect  of  successful  semi-automation  of  the  TCD  process 
but  recognize  that  we’ve  a  ways  to  go  before  we  succeed.  In  particular,  our  military  TCD  consultant,  Alex 
Razumovsky,  PhD,  has  verified  the  importance  of  this  approach  to  TCD  measurements  and  has  asked  us  to 
(a)  increase  the  depth  of  field  of  search  for  fastest  blood  flow,  to  (b)  improve  the  signal-to-noise  ratio  of  our 
system,  to  (c)  explore  use  of  an  intermediate  variable  -  a  time-varying  map  of  the  MCA  based  upon  creation  of 
an  ‘iso-surface’  of  Doppler-derived  blood-flow  speeds  -  a  ‘Doppler  map’,  discussed  below  -  to  map  the  blood¬ 
vessel  structure  sufficiently  to  visualize  the  MCA  in  vasospasm.  This  last  issue  is  particularly  important 
because  some  patients  after  blast  show  large  blood  flow  speeds  in  the  MCA  not  because  they  are  in 
vasospasm  but  because  their  hematocrit  is  sufficiently  low  that  the  resulting  reduced  effective  viscosity  of  the 
blood  allows  for  higher  than  normal  blood  flow  speeds. 

We  have  one  paper  in  print,  and  now  have  submitted  two  more  papers  that  describe  and  exercise  our 
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flow  envelope  calculation. 


Specific  Aim  #2:  To  modify  as  necessary  the  TCD  holder  for  use  on  the  ICA,  and,  the  algorithms 
required  to  localize  as  well  as  extract  the  highest  blood  flow  velocities  in  the  MCA  and  ICA,  in  order  to 
guarantee  that  the  aTCD  device  can  detect  the  highest  levels  of  vasospasm.  (Q3-Q12  out  of  12  quarters) 

Task  8:  Modify  as  necessary  the  brace  that  ordinarily  holds  the  transducer  to  the  temple  (to  facilitate 
measurement  of  blood-flow  from  within  the  MCA)  so  that  it  will  also  automatically  extract  blood  flow  information 
from  the  ICA  of  patients.  (Q4-Q8)  PSI 

(8.1)  Our  third  version  of  the  headset  (reported  in  our  second  Annual  Report)  now  covers  both  the  MCA 
and  internal  carotid  arteries.  Patients  remain  uncomfortable  using  any  headset,  however.  We  have  now 
developed  a  new  headset,  which  is  intended  to  minimize  patient  discomfort  while  stabilizing  the  transducers. 
Below  are  the  computer  modeled  plans  for  the  headset.  The  two  transducers  are  clipped  in  place,  and  can  be 
adjusted  by  sliding  the  holder  and  securing  it  with  the  knob.  This  should  allow  the  headset  design  to  function  on 
most  head  shapes  and  sizes.  Disposable  foam  pads  would  line  the  headset  and  would  be  replaced  for  each 
patient.  All  the  plastic  used  in  the  headset  should  be  autoclave  safe,  such  that  the  entire  headset  (except  for 
the  foam  pads)  could  be  autoclaved. 
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Figure  6. 

(A)  Computer  generated  model  of  the  final  headset 
design.  Both  transducers  are  placed  on  sliders  such 
that  they  can  be  accurately  positioned  for  a  given 
patient’s  anatomy. 

(B)  Another  view  of  the  computer  generated  model  of 
the  headset,  showing  the  placement  of  the  foam  pads 
between  the  patient’s  head  and  the  plastic  components 
of  the  headset. 

(C)  The  actual  produced  prototype  of  the  headset, 
modelled  on  a  real  person.  The  transducers  would 
wrap  around  the  back  of  the  head  to  avoid  interfering 
with  a  patient’s  movement. 


Milestone  #4:  Deliver  device  that  holds  the  transducer  that  supports  the  measurement  of  blood  flow 
velocity  in  the  MCA  and  ipsilateral  ICA  for  its  deployment  in  Aim  #1.  Done.  We’ve  encountered  an  important 
practical  problem  for  this  study,  however,  namely  that  no  test  subjects  wish  to  use  the  headset,  given  the 
general  discomfort  of  our  patients,  all  of  whom  have  not  had  sufficient  subarachnoid  hemorrhage  to  require 
induction  of  a  therapeutic  coma.  The  headset  is  comfortable  enough,  as  tested  by  engineers.  The  vasospasm 
patients  are  profoundly  uncomfortable,  however,  and  not  obliged  to  wear  the  headset  per  our  existing  IRB 
protocol.  We  will  talk  to  the  program  manager  about  adding  a  new  Task  to  this  study  where  we  will  refine  the 
designs  of  this  headset  using  feedback  from  healthy  volunteers. 

Task  9:  Modify  as  necessary  the  software  within  the  aTCD  device  to  optimize  its  ability  to  calculate  the 
value  of  fastest  blood  flow  (the  spectral  envelope)  within  the  MCA  and  the  ICA  for  vasospasm  patients.  (Q4- 
Q10)  UW  and  PSI. 

First  recall  that  the  device  as  delivered  will  report  values  of  the  spectral  envelope.  Our  job  requires 
optimization  of  this  process  for  vasospasm  patients,  a  notoriously  difficult  cohort.  We  have  identified  two 
pathways  towards  development  of  optimized  means  of  calculating  the  value  of  fastest  blood  flow  in  the  MCA 
and  ICA,  referred  to  here  as  the  spectral  envelope:  a  statistical  analysis  and  a  mathematical  modeling  analysis. 
We  have  published  results  from  the  first  pathway  and  continue  to  refine  the  details  of  the  second  method.  We 
have  now  begun  the  process  of  initiating  data  collection  on  vasospasm  patients  at  Swedish  Medical  Center 
here  in  Seattle  using  a  certain  standard  TCD  device  that  allows  us  to  access  the  raw  flow  data.  We  have  had 
difficulties  recruiting  enough  patients  that  actually  enter  vasospasm,  so  this  expansion  will,  hopefully,  lead  to  a 
greater  vasospasm  sample  size  on  which  to  test  our  enveloping  algorithms.  We  hope  to  begin  collecting  data 
on  vasospasm  patients  at  Swedish  by  the  end  of  2014. 

(9. 1)  As  reported  earlier  we  have  applied  our  Double-Gaussian  spectral  flow  envelope  calculator 
(Marzban  et  al,  2013)  to  our  new  vasospasm  patients  with  success  thus  far.  We  have  completed  work  on 
another  means  of  calculating  the  spectral  envelope.  This  new  approach  uses  the  ‘mixture’  paradigm  of 
Marzban  et  al,  2013,  where  a  sum  of  functions  adds  together  to  describe  the  distribution  of  blood-flow  speeds  a 
given  instant  in  time,  from  which  the  point  of  fastest  blood  flow  is  derived  for  that  instant.  The  earlier  published 
work  assumed  that  the  instantaneous  distribution  of  flow  was  best  modeled  by  a  Gaussian  function.  Here 
(Marzban  et  al,  2016)  we  relax  that  assumption,  making  use  of  two  alternatives  -  a  ‘skewed’  Gaussian  function 
and  an  entirely  non-Gaussian  function  (Figure  7-8).  The  skewed-Gaussian  function  has  two  additional 
parameters  -  the  shape  parameters  for  the  signal  and  the  background  components.  The  non-Gaussian 
function  uses  a  kernel  method  to  estimate  the  flow  velocity  distribution.  We  found  that,  when  it  worked,  the 
non-Gaussian  function  produced  moderately  superior  envelopes  than  the  skewed-Gaussian,  however  it  did  not 
work  for  all  patients.  We  found  that  the  skewed  Gaussian  function  produced  better  envelopes  than  the  regular 
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Gaussian  function  for  the  majority  of  patients  and  comparable  envelopes  (i.e.,  no  worse)  in  some  patients. 


GMM 


Skewed-GMM 


KMM 


Figure  7.  The  distribution  (more  accurately,  the  probability  density  function)  for  the  modeled  signal  (blue 
curve)  and  the  background  (red  curve),  as  determined  by  mixture  models  -  the  Gaussian  Mixture  Model  of 
Marzban  et  al,  2013  (GMM,  left  panel),  Skewed-GMM  (middle),  and  Kernel  Mixture  Model  (KMM,  right).  The 
blue  vertical  lines  denote  the  90th;  95th;  99th  percentiles  of  the  signal  flow  velocity,  and  the  black  vertical  line 
marks  the  maximum  FV  according  to  the  industry  standard  Modified  Geometric  Method. 


time  time  time 

Figure  8.  Sample  spectrograms  from  our  recently  published  paper  (Marzaban  et  al.  2016),  showing  the  fit  of 
envelopes  at  the  90th,  95th,  and  99th  percentile  of  signal  flow  velocity  for  each  of  the  Gaussian  Mixture  Model  of 
(GMM,  left),  Skewed-GMM  (middle),  and  Kernel  Mixture  Model  (KMM,  non-Gaussian,  right). _ 


(9.2)  We  have  submitted  our  paper  (Morison  et  al,  2015)  on  yet  another  means  of  calculating  optimal 
spectral  envelopes,  one  based  upon  a  first  principles  method  described  in  our  previous  report  with  final  results 
discussed  in  Figure  9-14,  below. 
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Figure  9  (left)  Comparison  of  measured 
(horizontal  axis)  versus  predicted  (vertical 
axis)  flow  envelope  for  each  of  the  five 
patients  we  studied  in  detail  in  Morison  et  al, 
(submitted).  There  is  decent  to  good  fit  on 
average,  with,  however,  a  tendency  of  our 
model  to  over-predict  the  measurements, 
(above)  Across  many  patients,  the 
regression  coefficient  averages  0.7  with  a  full 
distribution  noted  above. 


Here  we  use  blood  flow  modeling  because  it  helps  us  understand  the  underlying  flow  dynamics  and 
provides  diagnostic  value.  Spectrograms  can  often  be  corrupted  by  errors,  such  as  limits  associated  with  a 
finite  pulse  repetition  frequency  (setting  a  cap  on  the  largest  measurable  velocity),  Doppler  signal  from  multiple 
portions  of  the  vessel,  and  others.  Therefore,  if  we  can  apply  a  model  to  these  data  we  can  improve  the 
estimates  of  maximum  flow.  Figure  10  (below)  shows  the  major  components  of  the  method  for  modeling  flow. 
The  basic  idea  is  that  a  recursive  estimation  function  cycles  through  the  actual  collected  spectrogram  data  and 
the  predicted  spectrogram,  and  compares  the  two,  iterating  many  times  to  improve  the  agreement  between 
observed  and  predicted  data. 
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Figure  11.  A  sample  spectrogram 
used  as  input  to  the  first  principles 
model  of  Morison  et  al. 
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Figure  12.  The  output  flow  velocity,  where 
horizontal  arrow  length  is  proportional  to  the 
axial  velocity,  and  vertical  arrow  length  is 
proportional  to  the  radial  velocity.  The  artery  is 
shown  in  the  center  vertically,  and  the  x  axis 
represents  time. 


Below  is  an  example  of  an  observed  and  a  predicted  spectrogram,  using  the  model  we  have  developed. 
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Figure  13(A,B)-  An  observed  spectrogram  (left)  and  a 
predicted  spectrogram  (right). 

Figure  13(C).  The  correlation  between  our  data  and 
the  model.  Each  dot  represents  a  dataset  -  the  vertical 
coordinate  represents  the  correlation  between  the 
observed  data  and  our  methodology  and  the 
horizontal  coordinate  represents  the  correlation 
between  observed  data  and  the  null  methodology. 
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Figure  14.  Observed 
spectrograms  (left)  and  predicted 
spectrograms  (right)  for  which 
our  proposed  methodology  does 
best.  Each  pair  of  spectrograms 
represents  data  from  the  same 
patient. 


Milestone  #5:  Deliver  software  that  robustly  measures  the  values  of  fastest  blood  flow  measured  within 
the  MCA  and  ICA  for  testing  and  validation  in  Aim  #1.  UW.  As  described  above  we  have  developed  two 
plausible  pathways  to  replace  existing  commercial  algorithms  for  calculating  the  fastest  blood  flow  with  new 
ones  optimized  for  vasospasm,  where  the  signal  to  noise  ratio  is  low  and  the  blood  flow  speeds  are  high,  each 
testing  the  limits  of  the  existing  system.  In  future  work  we  will  seek  additional  means  of  calculating  this 
quantity  as  well  as  test  these  new  methods  on  our  existing  data  sets. 

Task  10:  Modify  as  necessary  the  TCD’s  controlling  software  so  as  to  optimize  its  ability  to  search  for 
and  locate  the  point  of  fastest  blood  flow  within  the  MCA  and  ICA  for  vasospasm  patients.  (Q4-Q1 0)  UW;  PSI 
(10.1)  In  year  2,  we  achieved  a  first  iteration  of  this  Task.  Our  studies  with  patients  made  clear  the  need 
to  have  the  ‘search  algorithm’  look  for  points  of  fastest  blood  flow  deeper  than  the  standard  60  mm  within  each 
of  the  middle  cerebral  artery  and  the  ICA.  We  have  made  that  change  (Figure  15,  reported  earlier)  and  reduced 
the  noise  in  the  collected  data. 
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Figure  15.  This  figure  shows  the  modified  user 
interface  that  guides  the  user  to  setup  the  automatic 
TCD  to  locate  the  point  of  fastest  blood  flow  in  the 
middle  cerebral  artery.  Here  we  photo-document 
(with  a  smart-phone  camera)  that  the  device  can 
search  up  to  70  mm  deep  within  the  cranium.  The 
user  specifies  the  search  depth  (60  or  70  mm)  and 
then  lets  the  machine  collect  its  data. 


Despite  the  fact  that  our  algorithm  can  search  over  a  depth  range  up  to  70  mm  within  the  cranium,  our 
search  and  locate  algorithm  thus  far  tends  to  lock  on  to  plausible  blood-flow  maxima  within  the  MCA  at 
shallower  depths  (40-50  mm)  as  compared  to  the  depth  of  vasospasm  identified  by  the  commercial  TCD  (50- 
65  mm),  which  facilitates  search  by  hand  for  the  global  maximum  in  blood  flow  speed  along  the  MCA.  In  Q9 
we  identified  a  potential  solution  to  this  issue  and  worked  with  our  consultant  Tom  Anderson  of  Glacier  River  in 
Q10  to  develop  the  necessary  software  changes.  In  Q1 1  we  updated  the  search  algorithm  to  allow  the  user  to 
specify  a  search  volume.  Here,  we  have  added  a  slider  bar  on  the  far  right  of  the  “Locate  Flow”  window  that 
allows  the  user  to  slide  a  search  box  spanning  10  mm  depth  across  the  accessible  volume  (below,  Figure  8a, b) 
We  can  access  volumes  between  40  mm  and  70  mm.  The  search  volume  box  can  be  moved  in  increments  of 
1  mm.  We  can  move  the  search  volume  box  for  both  the  MCA  and  ICA  transducers.  We  continue  to  advance 
another  portion  of  the  search  algorithm  design  to  improve  its  means  of  identifying  the  point  of  fastest  blood  flow 
within  a  given  volume. 
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Figure  16(A).  Screenshot  of  the  slider  (far  right)  and  search  volume  box  (center,  green  lines  bound  the  search 
box)  when  locating  flow  in  the  ICA. 


18 


Swedish  Hospital 


No  Patient 

Q 

Patient 


Step  1  -  Find  the  Temporal  Window 
Maneuver  the  transducer  in  the  general 
area  delineated  in  the  diagram,  and 
locate  the  position  of  maximum  signal 
strength  as  shown  by  the  indicator: 
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Step  2  -  Center  the  Flow 

Adjust  the  transducer  angle  to  target  the  vessel  and 
approximately  center  the  peak  flow  in  the  field  of  view: 
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Figure  16(B).  Screenshot  of  the  slider  (far  right)  and  search  volume  box  (center,  green  lines  bound  the  search 
box)  when  locating  flow  in  the  MCA. _ 


(10.2)  We  started  to  determine  what  it  will  take  to  access  an  internal  variable  within  the  aTCD  system 
that  facilitates  identification  of  the  point  of  fastest  blood  flow  within  the  cranium  in  the  region  of  the  MCA,  by 
producing  a  time-varying  map  of  an  iso-surface  of  Doppler-derived  blood-flow  speed  (Figure  17).  The  potential 
benefits  the  projects  are  three-fold.  (1)  We  expect  to  refine  the  process  of  identification  of  fastest  blood  flow  by 
analyzing  how  this  system  collects  and  processes  this  internal  imaging  data.  (2)  We  expect  to  produce  the 
ability  to  visualize  vasospasm  by  directly  measuring  the  diameter  of  the  blood  vessels  through  analysis  of 
these  Doppler  maps.  (3)  Production  of  these  maps  may  allow  assessment  of  vasospasm  in  the  case  of  low 
hematocrit  (which  can  produce  large  values  of  blood  flow  speed  without  vasospasmO  by  superimposition  of  the 
identified  point  of  fastest  blood  flow  with  the  associated  vascular  anatomy. 

This  will  require  substantial  work  ‘under  the  hood’  to  make  these  data  readily  available  for  both  off-line 

and  real-time  analysis.  We  have  started  the  process  of  setting  these  internal  tasks. _ 

Figure  17.  Samples  of  instantaneous  three- 
dimensional  maps  of  blood  flow  in  the  cranium. 
Here  we  display  in  one  map  a  three-dimensional  view 
of  the  internal  carotid  artery  (in  red)  and  where  it 
meets  the  middle  cerebral  artery  (MCA,  in  blue). 

Each  color  displays  a  single,  absolute  value  of  blood 
flow  speed  towards  or  away  from  the  transducer, 
thereby  outlining  the  blood  vessel.  The  other  three 
images  offer  three  different  perspectives  on  the  MCA 
(in  red)  and  the  anterior  communicating  artery  (blue), 
as  well  as  locating  the  position  of  the  point  of  fastest 
blood  flow  towards  the  transducer  (green  dot). _ 


TCD-angio  identified! 

(may  be  able  to  skip  CT-  or  MRI-angio 


we  can  map  blood 
vessel  structure  in  real 
time!  perhaps  can 
identify  vasospasm 
at  the  bedside 


(10.3)  We  improved  the  flow  detecting  algorithm  of  the  Presto  machine,  after  determining  that 
sometimes  the  peak  “power”  flow  site  detected  by  the  current  algorithm  was  not  the  desired  flow  site.  Ideally, 
the  algorithm  should  detect  the  peak  velocity  flow  site,  not  just  peak  power.  We  therefore  provided  two 
alternative  algorithms,  which  are  hybrids  of  power  and  velocity.  The  first  algorithm  searches  for  peak  power 
with  a  high  velocity  boost,  and  the  second  searches  for  peak  velocity  with  a  low  power  cutoff.  We  also  added  a 
Persistence  control  option,  which  is  accessible  in  all  three  of  the  imaging  modes  (Figure  18).  This  was  formerly 
hardcoded  in  the  Power  Imaging  mode. 
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Temporal  Window 

Maneuver  the 
transducer  in 
the  general 
area  delineated 
in  the  diagram, 
and  locate  the 
position  of  max 
signal  strength: 
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Locate  Flow  in  the  MCA  (Middle  Cerebral  Artery) 

Step  2  -  Center  the  Flow 

Ari|ust  the  transducer  angle  to  target  the  vessel  and 
approximately  center  the  peak  flow  in  the  field  of  view: 
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Figure  18. 

Implementation  of 
the  two  new  search 
algorithms  on  the 
Presto  machine. 
The  user  can  select 
between  the 
original  power 
imaging  mode  and 
the  two  new 
algorithms. 


The  power  with  high  velocity  boost  algorithm  is  computed  as  follows  (Figure  19,  below): 

if(  velocity  >  threshold  for  velocity  boost)  AND  (power  <  low  power  threshold)  AND 
(power  is  in  the  same  direction  as  the  search  direction) 

then  the  flow  estimate  is  power  multiplied  by  the  boost  factor 
otherwise 

the  flow  estimate  is  power  multiplied  by  the  derate  factor 
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Locate  Flow  in  the  MCA  (Middle  Cerebral  Artery) 

Step  2  -  Center  the  Flow 

Adfust  the  transducer  angle  to  target  the  vessel  and 
approximately  center  the  peak  flow  in  the  field  of  view: 
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Figure  19.  When 
the  Presto  is  using 
the  power  with  high 
velocity  boost 
search  mode, 
controls  are 
provided  specific  to 
that  algorithm. 


The  velocity  with  low  power  cutoff  algorithm  is  computed  as  follows  (Figure  20,  below): 

if  (power  >  low  power  cutoff)  AND  (power  is  in  the  same  direction  as  the  search  direction) 
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then  the  flow  estimate  is  velocity  multiplied  by  the  color  map  gain 
otherwise 
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Locate  Flow  in  the  MCA  (Middle  Cerebral  Artery) 

Step  2  -  Center  the  Flow 

Adjust  the  transducer  angle  to  target  the  vessel  and 
approximately  center  the  peak  flow  in  the  field  of  view: 
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Figure  20.  When 
the  Presto  is  using 
the  velocity  with  low 
power  cutoff  search 
mode,  controls  are 
provided  specific  to 
that  algorithm. 


Milestone  #6:  Deliver  software  that  optimally  identifies  the  point  of  fastest  blood  flow  in  the  MCA  and 
ICA  for  vasospasm  patients  for  testing  and  validation  for  its  analysis  in  Aim  #1.  UW 

Task  1 1 :  Perform  pilot  study  involving  32  patients  that  tests  in  a  blinded  fashion  the  efficacy  of  the 
aTCD  compared  with  a  standard  TCD  system.  (Q10-Q12).  UW 

(1)  We  present  these  results  in  Figures  4(A)  and  5  above,  showing  agreement  (within  +/1  12%) 
between  TCD  measurements  generated  via  our  research  TCD  device  and  a  clinical  TCD  device. 

Task  1 2:  Travel  to  a  military  medical  facility  under  the  sponsorship  of  Dr.  Armonda  to  determine  what 
changes,  if  any  need  to  be  made  to  the  FDA  approved  aTCD  system  for  optimal  deployment  to  military 
facilities.  Q10-12).  UW 

(1)  We  instead  facilitated  visits  by  Dr.  Razmovsky,  whose  input  is  evidence  in  the  headset  design  and 
in  our  pursuit  of  a  means  of  mapping  intra-cerebral  major  blood  vessels. 

Task  13:  Write  up  results  in  a  final  report  for  the  sponsor,  for  publication,  and  within  a  follow-on 
proposal  that  will  field  trial  an  aTCD  device  at  a  military  hospital  involving  patients  with  blast-induced  TCV. 
(Q10-Q12).  UW 

(1)  In  the  works. 

Milestone  #7\  Report  results  to  the  sponsor  and  begin  the  transition  of  this  technology  to  the  military. 
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KEY  RESEARCH  ACCOMPLISHMENTS 


•  We  have  created  an  improved  headset  for  use  in  deploying  simultaneous  MCA  and  ICA  monitoring. 

•  We  have  addressed  a  number  of  technical  obstacles  in  the  Presto  system  by  implementing  two  new 
methods  for  locating  the  peak  flow. 

•  We  have  determined  a  useful  enveloping  algorithm  that  tracks  well  with  clinical  TCD  data  and 
shows  improvement  in  the  Presto’s  peak  flow  detection  over  time. 

•  We  have  generated  a  useful  technique  for  modeling  blood  flow  based  on  actual  raw  spectral 
Doppler  data. 

•  In  total,  we  have  collected  119  data  sets  on  83  patients  over  the  3  years  of  patient  data  collection. 
We  believe  these  low  numbers  occurred  due  to  a  paucity  of  patients  to  consent  relative  to  historical 
averages,  and  because  only  one  third  of  the  available  patients  are  willing  to  participate  due  to  their 
extreme  discomfort. 
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CONCLUSIONS 


•  We  can  collect  high  quality  TCD  data  using  the  aTCD  system  from  vasospasm  patients. 

•  TCD  data  collected  on  the  Presto  tracks  well  with  standard  clinical  TCD  data  when  our  new 
enveloping  method  is  applied. 

•  We  have  generated  an  improved  TCD  headset  of  direct  use  for  our  aTCD  system  and  of  likely  use 
for  traditional  TCD  measurements. 
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Abstract 

Visually  evaluating  spectrograms  produced  by  Doppler  ultrasound 
has  proven  clinical  utility,  as  a  measure  of  blood  flow.  While  a  phys¬ 
iological  understanding  underlies  such  evaluations  an  explicit  model 
capable  of  approximating  any  given  spectrogram  has  not  been  avail¬ 
able.  Here,  we  propose  a  novel  physiological  model  in  which  a  simple 
circulatory  system  drives  pulsatile  flow  through  an  elastic  tube  and  a 
system  of  equations  relate  the  resulting  flow  velocity  field  to  their  cor¬ 
responding  simulated  spectrogram.  The  parameters  of  the  model  are 
then  optimized  to  match  measured  data  from  Transcranial  Doppler 
ultrasounds  to  produce  simulated  flow  patterns  and  thus  provide  a 
quantitative  physiological  interpretation  of  the  spectrogram. 

1  Introduction 

Modeling  blood  flow  has  been  a  topic  of  extensive  research  [1] ,  [2] ,  [3]  and  [4], 
It  is  important  not  only  for  the  purpose  of  understanding  the  underlying 
dynamics  [5]  and  [6],  but  also  for  its  diagnostic  value  [7]  and  [8]. 

Transcranial  Doppler  Ultrasound  (TDU)  provides  one  method  for  measur¬ 
ing  blood  flow  [9]  through  the  middle  cerebral  artery.  The  measurements  are 
generally  displayed  in  the  form  of  a  spectrogram.  A  spectrogram  is  a  measure 
of  the  intensity  of  acoustic  return  to  TDU  as  a  function  of  time  and  velocity 
1  Such  spectrograms,  however,  are  generally  corrupted  by  a  number  of  errors, 
including  limits  associated  with  a  finite  pulse  repetition  frequency  (which  sets 

1In  common  usage  the  term  spectrogram  refers  to  a  function  of  time  and  frequency;  but 
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a  cap  on  the  largest  measurable  velocity),  Doppler-signal  return  from  multi¬ 
ple  portions  of  the  blood  vessel  of  interest,  inadequate  acoustic  windows,  and 
other  sources  of  error.  A  basic  assumption  in  examining  Flow  Velocity  (FV) 
data  is  that  there  exists  an  underlying  spectrogram  which  adheres  to  certain 
physical  constraints.  For  example,  at  a  given  time,  the  brightness/amplitude 
of  a  spectrogram  is  expected  to  drop  to  zero  at  a  well-defined  FV  value  cor¬ 
responding  to  a  maximum  FVmax.  However,  observed  spectrograms  usually 
do  not  reflect  such  behavior.  Also,  observed  spectrograms  are  often  contami¬ 
nated  by  random  noise,  i.e.,  random  errors  which  cannot  be  attributed  to  any 
physical  constraints.  A  physics-  (or  physiologically)-based  model  is  required 
to  allow  one  to  attribute  physical  meaning  to  an  observed  spectrogram. 

The  approach  proposed  in  this  paper  involves  both  statistical  and  physi¬ 
cal  models.  Statistical  models  usually  involve  parametric,  mathematical  ex¬ 
pressions  whose  structure  is  often  dictated  by  the  convenience  of  estimating 
their  parameters.  Regression  models  are  canonical  examples  of  such  models 
(  [11]  and  [13]).  By  contrast,  physical  models  begin  with  physiologically- 
relevant  components  whose  mathematical  properties  are  well-understood, 
and  describe  the  overall  phenomenon  (e.g.,  blood  flow)  by  an  assembly  of 
such  components.  The  awkward  and  often  ill-conditioned  problem  of  fitting 
the  free  parameters  of  a  physical  model  to  data  is  referred  to  as  an  inverse 
problem.  2  Physiologically  based  blood  flow  models  of  the  type  examined 
here  have  been  considered  in  the  past,  but  the  model  parameters  have  not 
been  estimated  to  optimize  the  difference  between  observed  and  predicted 
spectrograms.  For  example  Mirsa  et  al.[8]  models  blood  flow  through  a  con¬ 
striction,  but  at  no  stage  in  that  work  is  a  spectrogram  used  to  estimate  the 
model  parameters.  Similarly,  Garcia  et  al.[10]  uses  a  physiological  model  to 
infer  a  velocity  field  from  Doppler  ultrasound  but  use  a  variable  focus  loca¬ 
tion  in  a  relatively  large  area  rather  than  a  spectrogram.  Notable  examples 
of  paired  physiology  and  ultrasound  models  are  given  by  Rune  Aaslid  and 
David  H  Evans.  Aaslid  offers  [19]  a  sophisticated  but  proprietary  simulation 
and  teaching  tool  and  Evans  employs  a  model  of  Doppler  ultrasound  in  mul¬ 
tiple  physiological  scenarios.  So  precedents  exist  for  the  forward  model  used 
in  this  work,  however  those  models  are  not  fit  to  observed  data.  When  a 

in  the  present  case  there  is  a  one  to  one  correspondence  between  frequency  and  Doppler 
velocity.  Also  a  spectrogram  is  a  function  of  time  in  the  sense  that  it  is  recorded  over 
time. 

2  The  distinction  made  here  between  statistical  and  physical  models  is  fuzzy,  but  ex¬ 
treme  examples  are  [14], [13]  and  [15], [16],  respectively. 
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spectrogram  is  fully  employed  as  an  objective  function,  in  for  example  Fer¬ 
nando  et  al.  [11]  and  Marzban  et  al.[13],  the  focus  of  the  work  is  an  envelope 
not  a  2-D  flow  velocity  field.  One  novel  feature  of  the  proposed  method  is  the 
manner  in  which  a  full  physiologically-based  blood  flow  model  is  utilized,  and 
its  parameters  are  estimated  from  an  entire  spectrogram.  Part  of  the  estima¬ 
tion  method  is  a  simple  instance  of  Recursive  Bayesian  Estimation  (RBE)[17] 
which  involves  iteratively  making  observations,  and  revising  estimates  about 
the  system  from  which  those  observations  originate.  RBE  has  been  proven 
useful  in  a  wide  range  of  applications  [18],  including  the  modeling  of  blood 
flow  [20]  and  [21]. 


Figure  1:  In  the  proposed  methodology  the  forward  model  consists  of  several  components,  each  with 
parameters  which  are  to  be  estimated  from  observed  spectrograms.  The  parameters  are  initialized,  a 
recursive  Bayesian  estimator  picks  a  time  series  through  the  output  of  the  forward  model  which  best 
matches  the  data.  Then  the  predicted  spectrogram  is  compared  with  the  observed  spectrogram  and  this 
whole  process  iterates  with  the  goal  of  improving  the  agreement  between  observed  and  predicted. 


The  flow  chart  in  Figure  1  highlights  the  major  components  of  the  method¬ 
ology.  As  shown  the  forward  model  is  comprised  of  three  sub-models.  First, 
a  Windkessel  model  3  in  the  style  of  [22]  produces  a  one  dimensional  time 
series  of  FVmax.  Second,  solutions  to  Navier-Stokes  equations  in  an  elastic 
tube  developed  by  [1]  and  [2]  are  set  to  match  the  FVmax  boundary  conditions 

3i.e.  an  electric  circuit  analogy  to  the  circulatory  system 
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and  return  a  2-D  flow  field  time  series.  Third,  an  equation  is  proposed  to 
calculate  the  spectrogram  resulting  from  measurement  of  the  2-D  flow  field 
with  TDU.  With  arbitrary  free  parameter  values  the  forward  model  outputs 
realistic  simulated  spectrograms,  but  the  output  does  not  resemble  any  par¬ 
ticular  patient.  The  purpose  of  the  methodology  developed  here  is  to  invert 
the  forward  model  and  find  input  parameters  corresponding  to  the  observed 
spectrogram  of  any  patient.  Most  of  the  free  parameters  are  constant  for  a 
given  patient.  For  example  neither  the  viscosity  of  blood  (p)  nor  the  angle 
of  the  transducer  relative  to  the  artery  (<f>)  are  assumed  to  change  within  a 
data  set.  The  phase  of  a  patient’s  cardiac  cycle  ( 9 )  is  expected  to  constantly 
advance  in  time  albeit  irregularly.  Given  values  for  all  other  free  parame¬ 
ters  the  forward  model  specifies  a  spectrogram  for  any  value  of  9.  The  RBE 
combines  an  observed  spectrogram  recorded  over  time  with  a  predicted  spec¬ 
trogram  modeled  over  the  domain  of  9  in  order  to  output  a  predicted  time 
series  of  9.  With  the  predicted  time  series  9  the  predicted  spectrogram  is 
compared  directly  to  the  observed  spectrogram.  This  comparison  feeds  back 
into  refined  values  for  the  free  parameters.  The  whole  process  repeats,  seek¬ 
ing  a  forward  model  that  best  matches  the  patients  data.  Seeking  to  invert 
our  model  with  this  methodology  many  iterations  of  the  loop  are  performed 
with  ad  hoc  shortcuts  and/or  in  parallel. 

2  Data 

Ninety-six  recorded  spectrograms  are  examined,  73  of  which  are  collected 
with  a  Spencer  Technologies  PMD  100  TDU  device  (Spencer  Technologies, 
Inc.,  Seattle  WA)  at  the  following  hospitals:  Harborview  Medical  Center 
(Seattle,  Washington),  Columbia  University  Medical  College  (New  York  City, 
New  York)  and  the  University  of  Texas  Medical  School  at  Houston,  Texas. 
The  remaining  23  spectrograms  are  collected  from  patients  at  Harborview 
Medical  Center  with  the  Presto  TDU  device  currently  under  development 
(PhysioSonics,  Inc.,  Redmond  WA).  Of  the  23  spectrograms  collected  with 
the  Presto  device  the  data  sets  are  anonymized  (as  are  all  of  the  data  sets) 
but  a  few  appear  to  be  drawn  from  repeated  patients,  suggesting  fewer  than 
23  separate  patients. 

All  of  the  spectrogram  data  analyzed  are  archived  at  a  sample  rate  of  40 
Hz  and  with  128  velocity  bins  from  a  bin-center  at  0  cm/s  to  a  bin-center 
at  308  cm/s.  All  patients  experienced  traumatic  brain  injury  (TBI),  with 
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the  majority  (59  out  of  the  73  Spencer  data  sets)  closed  head  TBI  and  the 
remaining  experiencing  penetrating  injuries.  Additionally  five  of  the  patients 
also  exhibited  vasospasm. 

The  spectrograms  for  the  96  patients  are  complex  and  varied  across  pa¬ 
tients.  Although  the  spectrogram  for  almost  every  patient  appears  to  be 
consistent  with  what  one  might  expect  from  an  ideal  spectrogram,  no  simple 
ideal  spectrogram  emerges  as  a  viable  candidate  when  multiple  patients  are 
examined.  In  fact,  the  spectrograms  for  9  of  the  patients  data  collected  with 
the  Spencer  device  are  sufficiently  complex  to  have  failed  a  commonly-used 
envelope-detection  algorithm,  referred  to  as  the  modified  geometric  algorithm 
[11].  Figure  2  shows  the  spectrogram  as  measured  by  a  Spencer  TDU  de- 


time  <s) 


Figure  2:  The  observed  spectrogram  for  one  patient  illustrating  many  relevant  phenomena.  FVmax  is 
easy  to  pick  out  as  a  function  of  time  demarcating  the  acoustic  clutter  (i.e.  back  ground  static)  above 
from  the  signal  of  interest  below.  Aliasing  is  visible  as  periodic  bright  spots  near  the  top  of  the  velocity 
axis  lining  up  with  systole.  A  variable  brightness  or  shading  is  clear  within  the  signal  of  interest,  in  this 
case  bright  near  FVmax  and  dim  near  velocity  equals  zero.  Also  the  ”  blip”  right  after  two  seconds  is  an 
example  of  non-additive  noise. 


vice.  Note  the  varying  shading  in  the  “body”  (i.e.,  the  signal  of  interest  in 
the  lower  portion)  of  the  spectrogram,  in  contrast  to  the  relatively  constant 
background  static.  It  is  important  to  note  that  the  amplitudes  shown  in 
this  figure  are  in  the  dimensionless  units  of  nepers  (Np)  or  the  natural  log 
of  the  archived  amplitude.  Indeed,  all  of  the  following  analysis  is  performed 
on  the  logarithm  of  the  observed  spectrogram  amplitude.  The  logarithmic 
transformation  is  motivated  by  the  empirical  result  that  the  histogram  of  the 
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Figure  3:  The  circuit  diagram  corresponding  to  the  Windkessel  model  used 


logarithm  of  the  amplitudes  is  bell- shaped,  where  as  without  the  logarithm 
the  histogram  is  heavy-tailed. 

The  data  sets  have  been  limited  to  three  second  excerpts.  Three  seconds 
was  chosen  for  the  reason  that,  for  most  patients,  it  is  long  enough  to  show 
multiple  cardiac  cycles  yet  short  enough  to  display  the  details  of  each  cycle. 
Sampled  at  40Hz  this  leads  to  120  time  steps  considered  from  each  patient. 

3  Methods 

3.1  Windkessel  Model 

Driven  with  a  four  level  square  wave,  the  circuit  illustrated  in  Figure  3  can 
produce  realistic  time  series  of  FVmax.  The  use  of  a  four  level  square  wave 
is  a  computational  expedient,  rather  than  a  physically  based  assumption.  A 
more  physically  reasonable  Windkessel  model  of  the  heart,  with  it’s  valves 
and  nervous  system,  would  require  nonlinear  circuit  elements.  When  at¬ 
tempting  to  invert  the  forward  model  such  a  nonlinear  circuit  would  become 
a  computational  bottle  neck,  while  the  circuit  in  Figure  3  may  be  simulated 
at  negligible  computational  cost.  Free  parameters  associated  with  the  various 
components  of  the  Windkessel  model  are  tabulated  in  Table  1.  As  described 
in  text,  most  of  the  free  parameters  are  constant,  while  phase  in  cardiac  cycle 
(i 9 )  changes  with  time.  Table  1  is  written  to  present  the  true  count  of  free 
parameters.  In  practice  the  domains  of  the  free  parameters  are  constrained 
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in  an  ad  hoc  manner,  for  example  a  linear  combination  of  r  and  the  V’s  are 
adjusted  twice  per  iteration  of  the  outer  loop  in  1.  So  there  are  ten  constant, 
plus  one  dynamic,  free  parameters  but  in  practice  their  domains  are  inter¬ 
related.  An  example  of  the  Windkessel  model  output  is  shown  in  Figure  4. 


parameter 

symbol 

domain 

units  (abstract) 

phase  in  cardiac  cycle 

e 

0  to  2n 

rad 

heart  rate 

CO 

0.67T  to  67T 

i 

second 

resistance 

r 

0  to  1 

pressure 

flow 

capacitance 

c 

0  to  10 

flow-time 

pressure 

inductance 

l 

0  to  0.2 

pressure-time 

flow 

dwell  a 

A 

0.05  to  0.25 

time 

dwell  b 

A  b 

0.05  to  0.25 

time 

dwell  c 

A  c 

0.00  to  0.20 

time 

amplitude  a 

Va 

0.0  to  1. 

pressure 

amplitude  b 

vb 

0.0  to  1. 

pressure 

amplitude  c 

Vc 

1 

pressure 

Table  1:  The  dynamic  and  fixed  parameters  introduced  in  the  Circulatory  System  Calculation  component 
of  the  model.  The  physical  units  are  general  because  the  Windkessel  model  is  an  analogy,  ’’flow”  will  be 
scaled  to  an  FVmax  time  series  in  cm/s  before  it  is  input  to  the  flow  field  calculation. 


3.2  Flow  Field  Calculation 

The  time  series  FVmax  constitutes  the  boundary  condition  for  a  fluid  dy¬ 
namic  model  of  blood  flow  in  the  middle  cerebral  artery.  In  order  to  use 
solutions  developed  by  [1]  and  [2]  the  middle  cerebral  artery  is  treated  as 
an  elastic  tube,  and  blood  is  assumed  to  satisfy  the  incompressible  Navier- 
Stokes  equations.  The  tube  is  assumed  to  be  straight;  this  assumption  is 
reasonable  because  although  the  middle  cerebral  artery  is  convoluted,  the 
typical  radius  of  curvature  along  the  artery  is  larger  than  the  radius  of  the 
artery  itself.  The  vessel  wall  is  assumed  to  be  thin,  allowing  one  to  treat  it 
as  a  Hookean  membrane[23],  i.e.,  the  2-dimensional  generalization  of  an  ideal 
spring.  We  seek  solutions  for  radial  and  axial  velocities  of  the  form  v  (r,  t) 
and  u  (r,  t),  where  r  is  the  distance  from  the  axis.  As  such,  the  solutions  are 
rotationally  symmetric,  and  there  is  no  angular  velocity.  As  is  often,  done 
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Figure  4:  An  example  of  FVmax  produced  by  the  circuit  shown  in  Figure  3 


approximations  are  made  in  order  that  solutions  may  be  obtained  through 
a  Fourier  decomposition  of  the  forcing.  In  the  end,  two  degenerate  modes 
of  wave  propagation  are  found,  for  each  harmonic  component.  In  most  prior 
work  an  appeal  to  experiment  has  been  made  at  this  point  to  determine  the 
mode  of  wave  propagation.  In  our  case  the  relative  phase  and  amplitude 
between  the  modes  becomes  a  free  parameter.  The  forward  model  inversion 
embodied  in  the  flow  chart  of  Figure  1  becomes  an  effort  to  determine  the 
mode  of  wave  propagation  for  each  patient.  The  mechanical  properties  of 
blood  and  vessel  walls  are  the  subject  of  many  studies  [27]  [28]  [29]  [30]  [31] 
[32]  [33]  which  provide  useful  bounds  for  the  free  parameters  (Table  2)  within 
which  to  optimize  their  estimated  value  from  data.  The  output  of  model  of 
blood  flow  in  the  middle  cerebral  artery  is  shown  in  Fig.  5.  Typically,  the 
radial  velocity  magnitude  is  much  less  than  of  the  axial  velocity  magnitude; 
here,  they  have  been  scaled  to  achieve  visual  clarity.  It  may  be  tempting  to 
ignore  the  lower  amplitude  radial  component  in  favor  of  the  higher  amplitude 
axial  component;  however,  because  the  observed  spectrograms  span  the  full 
domain  of  velocity,  including  near-zero  velocity,  the  scale  disparity  between 
the  two  components  does  not  justify  discarding  the  lower  amplitude  veloci¬ 
ties.  The  mathematical  details  of  the  flow  model  developed  here  are  similar 
to  that  developed  in  [23]  4. 

4except  that  Equations  17.134  and  17.135  seem  to  be  in  error  and  the  reader  will  need 
to  derive  the  correct  equations  from  17.128,  17.129,  and  17.131 
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parameter 


symbol  domain 


units 


density  of  blood 
viscosity  of  blood 
density  of  vessel  wall 
equilibrium  radius  of  vessel 
thickness  of  vessel  wall 
Poisson’s  ratio  of  vessel  wall 
Young’s  modulus  of  vessel  wall 
relative  amplitude  of  ”  +”  mode 
relative  amplitude  of  mode 

Table  2:  The  fixed  parameters  introduced  in 


p 

1.06 

9 

cm 3 

p 

3  to  4 

centipoise 

Pw 

0.98 

9 

cm 3 

a0 

0.15  to  0.25 

cm 

h 

5  to  20 

o 

O 

Vw 

0.3  to  0.7 

dimensionless 

E 

40  to  100 

9 

cm-  s 

*+ 

unit  normal 

complex  plane 

Y_ 

unit  normal 

complex  plane 

3  Flow 

Field  Calculation  component  of  the  model. 

time  (s) 


Figure  5:  The  output  flow  velocity  is  shown.  The  horizontal  arrow  length  is  proportional  to  the  axial 
velocity  and  the  vertical  arrow  length  to  the  radial.  The  arrow  base  positions  increment  horizontally  with 
time,  and  vertically  with  radial  distance.  To  make  the  graphic  more  intuitively  clear  the  artery  center  is 
shown  in  the  center  vertically  and  reflected  axial  velocities  and  radial  positions  are  plotted. 
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3.3  Spectrogram  Calculation 

The  preceding  steps  in  the  proposed  model  produce  a  2-D  flow  velocity  field. 
TDU,  however,  produces  a  spectrogram,  i.e.,  a  time  series  of  the  histogram 
of  the  volume  corresponding  to  a  Doppler  velocity  v  where  volume  is  as¬ 
sumed  to  be  proportional  to  the  total  number  of  scatterers  [24].  In  order 
to  compare  the  simulated  blood  flow  velocity  field  to  the  observed  Doppler 
ultrasound,  therefore,  a  spectrogram  must  be  derived  from  the  modeled  ve¬ 
locity  field.  The  TDU  considered  here  interrogates  only  the  flow  velocity 


Figure  6:  An  overview  of  the  off-bore  angle  <f>  and  the  polar  angle  <f>.  This  geometry  determines  the 
Doppler  velocity  from  the  transducer’s  perspective  on  the  axisymmetric  flow  field. 


field  component  parallel  to  the  transducer  axis.  Though  the  simulated  blood 
flow  is  axisymmetric  the  transducer  is  not  assumed  to  be  aligned  with  the 
middle  cerebral  artery.  Introducing  an  off-bore  angle  <f>,  the  speed  toward 
the  transducer  v'  is  given  by 

v'  (r,  0,  <f>)  =  vaxiai  (r)  cos  <f>  +  vradial  (r)  sin  $  cos  0  ,  (1) 

with  0  the  angle  in  cross  section  of  the  artery  (see  Figure  6).  Said  differently, 
0  is  the  polar  angle  around  the  symmetry  axis  of  the  artery  and  $  is  the 
angular  misalignment  between  the  two  symmetry  axes,  that  of  the  ultrasound 
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focus,  and  that  of  the  artery.  The  location  of  the  TDU  focus  is  assumed  to 
be  offset  from  the  artery  center  and  variable  in  width.  A  finite  focus  volume 
enters  into  the  calculation  of  a  spectrogram  simply  as  a  location  dependent 
weighting  function  /  (r,  0).  Values  of  v'  from  equation  1  are  shown  in  Figure 
7.  The  flow  velocity  field  is  the  same  as  is  shown  in  Figure  5  but  in  that 
graphic  the  flow  is  displayed  over  time  and  a  diameter  of  the  artery.  This 
graphic  illustrates  the  effect  of  Tand  0  with  a  snapshot  of  v'  throughout  a 
cross  section  of  the  artery.  An  ideal  spectrogram  is  a  measure  of  how  much 


cm/s 


left  right  (cm) 


Figure  7:  The  Doppler  velocity  component, vf ,  is  shown  as  a  scaled  color  in  the  symmetry  plane  of 
the  artery.  As  shown  r  and  <j>  would  be  polar  coordinates  in  the  Cartesian  coordinate  axes,  with  </>  =  0 
illustrated  as  being  the  right.  For  this  figure  a  relatively  large  off-bore  angle  of  <E>  =  ^-rad  =  80°  is  chosen 
to  illustrate  the  (f)  dependence  of  equation  1.  In  cases  of  <E>  closer  to  0  the  axial  contribution  to  the  Doppler 
velocity  dominates  and  the  second  term  in  equation  1  is  not  apparent. 


volume  within  the  focus  area  is  moving  at  a  velocity  v  In  other  words,  given 
time  t  and  velocity  v  the  spectrogram  is  a  measure  of  the  infinitesimal  volume 
for  which  v'  =  v  as  a  function  of  time.  An  ideal  spectrogram  s^eai  ( v )  can  be 
defined  as 


s  (v)  =  lim  — Area  in  r,  0  such  that  \v'  (r,  0)  —  v\  <  e  (2) 

where  the  area  is  understood  to  be  either  inside  the  ultrasound  focus  or  to  be 
weighted  by  some  function  defining  the  focus  /  (r,  0).  The  ideal  spectrogram 
is  the  infinitesimal  area  for  which  the  Doppler  component  of  the  flow  velocity 
field  inside  the  ultrasound  focus  is  equal  to  the  ideal  spectrogram  velocity. 
It  is  relatively  straightforward  to  show  that  for  a  broad  class  of  v’ ,  the  limit 
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in  equation  2  can  be  written  as 


s(u) 


|  Vn'| 


(3) 


where  l  is  any  parameterization  of  contours  for  which  v' =  v.  The 
relationship  between  an  ideal  spectrogram  and  a  flow  velocity  field  has  been 
studied  by  Evans  [4],  but  not  in  integral/closed  form  as  in  equation  3.  The 
focus  /  (r,  4>)  is  defined  in  the  symmetry  plane  of  the  artery  and  specifically 
to  be  a  2-D  Gaussian  inside  of  the  artery  and  zero  outside. 


f  (r,  <t>) 


(r  cos  (j>  —  offsetx)^+(r  sin  <j>  —  offset^ 

e  2width^  /  (27rwidth2 ) 

0 


,r  <a0 
,r  >  a0 


(4) 


This  2-D  Gaussian  introduces  three  free  parameters:  offset.,,  an  offset  in  the 
symmetry  plane  of  the  artery  parallel  to  4>  =  0,  offset^  an  offset  in  the  sym¬ 
metry  plane  of  the  artery  parallel  to  0  =  y ,  and  width  a  parameter  governing 
the  width  of  the  focus.  When  width  is  large  i.e.  near  2ao,  the  offset  param¬ 
eters  have  little  effect  on  the  predicted  spectrogram,  reverting  to  the  case  of 
uniform  insonification.  In  practice  the  bloods  Doppler  velocity  v'  is  approxi¬ 
mated  with  a  piecewise  planar  triangular  surface  and  the  integral  in  equation 
3  is  approximated  with  a  sum.  This  is  one  of  the  more  computationally  ex¬ 
pensive  steps  in  the  methodology.  To  avoid  performing  this  calculation  twice 
per  iteration  of  the  outer  optimization  loop  in  the  flow  chart  Figure  1  the 
spectrogram  due  to  the  tissue  motion  is  calculated  with  a  quick  approxima¬ 
tion.  Additionally,  the  observed  spectrogram  is  assumed  to  differ  from  the 
ideal  spectrogram  by  the  addition  of  two  types  of  noise  and  by  the  unknown 
acoustic  reflectivity  of  blood  and  tissue. 


S  (l7,  9)  =  BbS  ( Vaxiai  ( 9 )  ,  Vradial  ( 9 ))  +  BtS  (f]  ( 9 ) ,  f  (6>))  +  B0  (5) 

Z(v,t)  =  S(v,9(t))  +  N(v,t)  (6) 

Here  S  is  the  idealized  spectrogram  output  by  the  forward  model  over  the 
domain  of  spectrogram  velocity  v  and  phase  in  the  cardiac  cycle  9.  Z  is  the 
observed  spectrogram  data  which  differs  from  S  by  zero  mean  one  dimen¬ 
sional  Gaussian  random  variable  (i.e.  noise)  N  and  an  unknown  relationship 
between  the  domain  variables  9  and  time.  The  free  parameters  B0,  Bb  and 
Bt  are  the  acoustic  brightness  of  the  background  clutter,  blood  and  tissue 
respectively. 
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parameter 


symbol 


domain 


units 


off-bore  angle  of  transducer 

$ 

tan-1 

radians 

vessel  wall  Doppler  spread 

0  to  3 

j_ 

dv 

focus  width 

focus.width 

0.1  to  2. 

J_ 

o0 

focus  offset  x 

focus  .offset  _x 

N 

o0 

focus  offset  y 

focus  .offset  _y 

N 

a,0 

reflectivity  of  blood 

Bb 

— oo  to  oo 

Np 

reflectivity  of  tissue 

Bt 

— oo  to  oo 

Np 

background  clutter 

Bo 

— oo  to  oo 

Np 

Table  3:  parameters  introduced  in  the  Spectrogram  Calculation  component  of  the  model.  The  abbrevi¬ 
ation  Np  stands  for  neper  of  the  spectrogram  amplitude.  A  neper  is  to  the  natural  log  as  a  decibel  is  to 
a  base  ten  logarithm. 


3.4  Recursive  Bayesian  Estimation  (RBE) 

The  forward  model  outputs  an  ideal  spectrogram  S  (v,6),  an  expected  am¬ 
plitude  as  a  function  of  velocity  and  phase  in  the  cardiac  cycle.  Of  course  the 
observed  spectrogram  z(v,t)  is  not  naturally  delineated  into  cardiac  cycles. 
Within  a  recursive  Bayesian  framework  the  observed  spectrogram  z(v,t )  is 
combined  with  S  (v,  9)  to  produce  a  posterior  probability  density  function 
(PDF)  over  6  which  evolves  in  time.  The  location  of  the  maximum  value  of 
this  posterior  PDF  is  taken  as  an  estimate  of  the  time  series  9  (t).  In  sum- 
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Figure  8:  On  the  left  is  an  example  of  the  likelihood  of  the  data,  from  one  patient,  conditioned  on  6  as  a 
function  of  time.  On  the  right  is  the  cumulative  probability  density  function  over  0,  as  a  function  of  time. 
For  the  earliest  time  steps,  the  cumulative  PDF  is  ambiguous  with  regards  to  the  value  of  0,  but  soon  a 
well  defined  mode  in  9  develops.  The  mode  advances  in  time  roughly  in  agreement  with  «  u. 

mary,  the  forward  model  produces  an  ideal  spectrogram  over  9.  The  RBE 
compares  each  time  step  of  the  observed  spectrogram  to  the  ideal  spectro- 
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gram  at  each  each  value  of  9.  Then  the  RBE  uses  the  information  that  9 
advances  in  time  at  approximately  the  heart  rate  (u)  to  predict  a  time  series 
of  9.  RBE’s  have  a  long  history  of  use  in  signal  processing  and  the  details  of 
any  given  implementation  may  vary  dramatically.  Stylistically  the  RBE  used 
here  is  similar  to  one  used  in  Nodestar  [18]  in  a  few  ways.  The  PDF  of  the 
state  variable  (9  in  our  case)  is  stored  as  a  relatively  fine  grained  histogram. 
The  presumed  relationship  between  the  state  variable  and  a  measurements 
of  the  system  is  not  linearly  approximated.  Also,  discretized  numerical  inte¬ 
gration  is  used  at  each  time  step  to  marginalize  over  previous  values  of  the 
state  variable.  These  attributes  are  all  in  contrast  to  the  most  typical  RBE, 
a  Kalman  filter.  Our  RBE  is  formulated  as  it  is  to  accommodate  a  relatively 
complex  relationship  between  9  and  S  (v,  9). 

3.5  Parameter  Refinement  and  Feature  Extraction 

Each  of  the  components  of  the  proposed  forward  model  introduces  a  number 
of  free  parameters  (Tables  1,2  and  3)  and  with  the  exception  of  9  all  are 
assumed  constant  throughout  a  patients  data  set.  With  a  complete  set  of 
estimated  free  parameters  including  9  (t)  the  forward  model  as  a  whole  is 
compared  to  the  data.  The  total  sum  squared  error  between  the  predicted 
spectrogram  and  the  observed  spectrogram  is  assumed  to  determine  the  joint 
probability  of  that  set  of  free  parameter  values  5.  In  addition  to  the  stated 
domain  over  which  each  free  parameter  will  be  optimized,  multiple  ad  hoc 
estimates  and  refinements  are  made,  graphically  implied  with  curved  green 
lines  radiating  from  ’’Data:  Observed  Spectrogram”  in  the  flowchart  Figure 
1.  For  example  the  acoustic  reflectivity  of  blood  (5&)  is  stated  to  be  opti¬ 
mized  between  —  oo  and  oo  because  each  successive  value  comes  from  a  linear 
estimate.  In  fact  B b  is  refined  twice  per  iteration  of  the  loop  in  1,  first  set 
equal  to  values  given  by  linear  estimation  depending  on  a  subset  of  the  for¬ 
ward  model  and  secondly  with  the  whole  set  of  free  parameter  trial  values. 
The  most  effective  ad  hoc  optimization  strategy  stems  from  the  observation 
that  the  value  of  9  (t)  given  by  the  RBE  is  quite  stable  or  consistent  for 
every  free  parameter  except  heart  rate  (cj).  To  exploit  this  stability  many 
iterations  of  the  loop  in  1  are  performed  in  parallel  all  sharing  a  compromise 

sThe  reader  may  object  that  the  variance  of  the  noise  itself  should  be  a  free  parameter 
subject  to  optimization,  and  while  this  would  not  directly  affect  the  optimum  set  of  free 
parameters  found  at  this  step  it’s  value  is  implicit  in  0  (t). 
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uj  and  0  (t).  6  Successful  optimization  to  a  global,  or  even  local,  best  set 
of  free  parameters  is  the  subject  of  future  work.  Our  current  optimization 
simply  searches  subsets  of  free  parameters  iterating  a  fixed  number  of  times 
and  returns  the  best  candidate  values. 

Armed  with  estimates  of  the  dynamic  parameters,  one  can  generate  a 
variety  of  useful  quantities,  such  as  a  spectrogram,  an  envelope,  mean  flow 
or  pulsatility  index.  An  example  of  a  predicted  spectrogram  is  shown  in 
Figure  9,  where  it  is  compared  with  the  corresponding  observed  spectrogram. 
The  fitted  spectrogram  closely  matches  the  observed  spectrogram  in  multiple 
ways.  For  example:  the  gross  features  such  at  systolic  and  diastolic  velocities 
are  similar,  the  overall  brightness  matches  and  both  have  a  subtle  dicrotic 
notch.  Also  at  the  time  of  systole  both  the  data  and  the  predicted  have  a 
slightly  less  bright  spectrogram  for  velocities  less  than  FVmax  as  compared 
to  times  of  diastole.  A  pulsatility  index  can  be  constructed  as  follows:  The 
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Figure  9:  An  example  of  an  observed  (left)  and  a  predicted  (right)  spectrogram. 


pulsatility  is  defined  [6]  over  a  cardiac  cycle  as  the  difference  between  the 
systolic  and  diastolic  flow  velocity,  divided  by  the  mean  flow  velocity,  where 
systolic  and  diastolic  are  synonymous  with  the  maximum  and  minimum  of 
flow  velocity.  As  any  set  of  fixed  parameters  defines  a  cardiac  cycle,  they 
also  define  a  pulsatility  index.  Similarly,  any  physical  quantity  that  can  be 
calculated  from  the  forward  model  is  thereby  given  by  the  fixed  parameters. 

6For  initial  iterations  of  the  outer  loop  in  Figure  1  the  Windkessel  model  is  used  with 
fast  approximations  to  the  flow  field  and  spectrogram  calculations  in  order  to  perform 
parallel  iterations  on  a  single  processor.  A  second  round  of  iterations  employ  the  full 
forward  model  with  the  flow  field  and  spectrogram  calculations  but  rather  than  repeatedly 
use  the  RBE  a  single  uj  and  6  time  series  is  shared  for  many  iterations.  Though  it  would 
be  relatively  easy  to  achieve  we  have  yet  to  take  advantage  of  hardware  parallelization, 
because  other  algorithmic  improvements  are  expected  to  yield  dramatic  speed  gains. 
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4  Results 


In  order  to  assess  the  success  of  the  above  methodology  it  is  applied  to 
all  of  the  data  sets.  First,  the  goodness  of  fit  is  assessed  quantitatively  by 
comparison  to  a  null  methodology.  Second,  the  output  is  judged  qualitatively, 
through  a  visual  comparison  of  the  predicted  to  the  observed  spectrogram. 

Though  there  are  notable  examples  of  simulated  spectrograms  [19]  and  [4] 
there  seem  to  be  no  instances  of  simulated  spectrograms  fit  to  data.  There 
are,  however,  well  established  algorithms  to  draw  the  maximum  flow  velocity 
time  series  from  an  observed  spectrogram.  Any  FVmax  fitting  algorithm 
may  be  considered  to  produce  a  fitted  spectrogram  with  the  following  simple 
heuristic.  Any  given  time  series  FVmax  segregates  the  observed  spectrogram 
into  pixels  which  are  due  to  blood  flow  (FV  <  FVmax)  and  pixels  which  are 
not  (FV  >  FVmax).  All  pixels  due  to  blood  flow  are  assumed  to  have  values 
drawn  from  a  normal  distribution  with  a  mean  value  of  nullsignai.  Similarly 
the  values  of  the  pixels  not  due  to  blood  flow  are  assumed  to  be  drawn  from  a 
normal  distribution  with  the  mean  value  nullciutter-  The  modified  geometric 
method  [11]  is  long  established,  and  is  currently  in  use  in  the  Spencer  device. 
In  the  null  method  the  modified  geometric  method  provides  FVmax  and  the 
sample  mean  of  the  pixels  above  FVmax  are  taken  as  an  estimate  of  null  Gutter 
and  the  sample  mean  of  the  pixels  below,  an  estimate  of  nullsignai  An  example 
of  this  null  methodology  spectrogram  is  shown  in  Figure  10.  The  pixel  by 


Figure  10:  An  example  of  a  spectrogram  predicted  with  the  null  methodology  (right)  along  with  the 
observed  spectrogram  (left) 


pixel  Pearson’s  correlation  between  predicted  and  observed  spectrogram  is 
calculated  for  each  patient.  This  is  done  both  for  spectrograms  predicted 
by  the  methodology  of  this  paper  and  also  for  spectrograms  predicted  by 
the  null  method.  Using  the  correlation  between  predicted  and  observed,  as 
a  simple  performance  measure,  we  compare  the  proposed  methodology  to 
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the  null  method  in  Figure  11.  This  within-data  set  correlation  shows  clearly 


correlation  between  model  and  data 


Figure  11:  Each  dot  represents  one  data  set.  The  correlation  between  the  observed  data  set  and 
the  spectrogram  predicted  by  the  methodology  of  this  paper  is  the  vertical  coordinate.  The  horizontal 
coordinate  is  the  correlation  between  the  observed  data  the  null  methodology 


that  the  proposed  methodology  more  closely  fits  the  data,  but  to  asses  the 
goodness  of  fit  we  need  to  account  for  the  number  of  free  parameters  in  the 
null  versus  the  proposed  methodologies. 

The  Akaie  information  criterion  (AIC)  is  a  commonly  used  measure  of 
goodness  of  fit,  [12]  intended  to  rank  the  generalizabiliy  and  performance  of 
statistical  models. 

AIC  =  2 Nparams  —  2  In  (P  (data\par corns))  (7) 

Where  the  number  of  free  parameters,  Nparams,  and  the  likelihood  of  observ¬ 
ing  the  data  given  the  optimized  parameter  values,  P  ( data\params )  give  a 
metric  which  is  predictive  of  a  model’s  cross  validation  performance.  The 
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smaller  a  model’s  AIC  the  better.  The  null  method  gives  an  independent 
value  of  FVmax  for  every  time  step  in  the  data,  giving  120  free  parameters. 
It  also  assigns  a  mean  spectrogram  value  for  pixels  above  and  below  FVmax, 
giving  two  more  free  parameters.  The  proposed  methodology  assigns  a  dif¬ 
ferent  phase  in  the  cardiac  cycle  for  each  time  step  in  the  data,  giving  120 
free  parameters.  The  Windkessel  model  introduces  ten  more  parameters,  the 
flow  field  calculation  nine  and  the  spectrogram  calculation  eight.  So  for  the 
record  length  considered  the  proposed  methodology  has  a  total  of  147  free 
parameters  per  data  set  and  the  null  method  122.  It  is  clear  from  figures 


Figure  12:  Each  dot  shows  the  AIC  for  a  patients  data  set,  the  null  method  horizontally  and  proposed 
method  vertically.  With  one  exception  every  dot  is  below  the  diagonal  indicating  the  proposed  method¬ 
ology  is  better.  The  data  set  corresponding  to  the  lone  point  above  the  diagonal  appears  to  consist  solely 
of  static  and  one  ’’blip”. 


11  and  12  that  the  proposed  methodology  improves  upon  the  null  method. 
However  the  null  method  is  a  minimal  standard.  The  true  significance  of  this 
work  is  better  demonstrated  via  a  qualitative  result. 
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A  ready  application  of  the  methodology  is  in  blind  signal  processing  of 
observed  spectrograms.  Since  there  is  no  ground  truth  reference  signal  the 
reader  may  compare  the  original  observed  spectrogram  in  the  left  column  of 
figures  13,  14,  15  and  16  to  the  predicted  on  the  right,  and  judge  subjectively 
how  well  the  signal  in  the  data  is  captured  by  the  predicted  spectrogram.  To 
convey  some  of  the  variety  present  in  the  96  data  sets  present  the  five  worst 
and  the  five  best  are  shown.  The  terms  ’’best”  and  ’’worst”  are  defined  by  the 
relative  performance  measure  difference  between  the  null  and  the  proposed 
methods.  This  is  done  for  both  performance  measures  leading  to  ten  data 
sets  for  consideration.  Notice  that  these  ’’worst”  and  ’’best”  data  sets  are 
highlighted  in  black  and  yellow  respectively  in  figures  11  and  12 

5  Summary  and  Discussion 

A  methodology  has  been  proposed,  modeling  blood  flow  through  an  artery 
and  the  measurement  of  it  with  ultrasound.  The  approach  involves  a  com¬ 
bination  of  physiologically-based  and  statistical  models.  For  a  given  patient, 
each  component  of  the  overall  model  involves  parameters  that  are  estimated 
from  observed  spectrograms.  The  resulting  models  may  then  be  used  to  pro¬ 
duce  a  variety  of  useful  quantities,  including  the  spectrogram,  the  envelope, 
and  the  flow  field. 

A  fundamental  signal  processing  limitation  in  use  of  TDU  leads  to  alias¬ 
ing,  namely  spurious  features  at  the  top  of  the  spectrogram  which  appear  to 
reflect  the  true  signal  of  interest  at  the  bottom  of  the  spectrogram.  Another 
artifact  that  signal  processing  is  designed  to  minimize  is  the  true  acoustic 
signal  due  to  reflections  from  the  vessel  wall,  which  tends  to  be  brighter  than 
the  signal  from  blood  flow.  Because  aliasing  and  the  acoustic  signal  from 
the  vessel  wall  are  explicitly  included  in  the  forward  model  they  are  trivially 
easy  to  remove  from  the  predicted  spectrogram.  Figure  17  shows  an  illustra¬ 
tion  of  this  utility.  The  bottom  left  figure  qualitatively  resembles  the  true 
spectrogram  (top);  the  bottom  right  figure  shows  the  spectrogram  with  the 
aliasing  and  vessel  wall  reflection  removed. 

Another  novel  feature  of  the  work  here  is  the  articulation  of  a  closed- 
form  relationship  between  a  general  two  or  three  dimensional  flow  field  and 
a  spectrogram,  namely  Equation  3.  This  expression  allows  one  to  compute 
a  spectrogram  from  any  flow  field,  regardless  of  how  it  is  generated.  The 
ability  to  relate  the  flow  field  to  a  spectrogram  is  important  for  any  attempt 
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Figure  13:  Observed  spectrogram  (left)  and  predicted  spectrogram  (right)  of  the  five  patients  for  which 
the  proposed  methodology  is  ” worst”.  In  this  case  we  mean  worst  in  the  sense  that  difference  between  the 
pixel  by  pixel  correlation  of  the  methodology  minus  the  null  method  is  minimal. 
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Figure  14:  Observed  spectrogram  (left)  and  predicted  spectrogram  (right)  of  the  five  patients  for  which 
the  proposed  methodology  does  ’’best” .  Where  we  mean  best  in  the  sense  that  the  correlation  of  the 
proposed  methodology  minus  correlation  of  the  null  method  is  maximal. 
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Figure  15:  Observed  spectrogram  (left)  and  predicted  spectrogram  (right)  of  the  five  patients  for  which 
the  proposed  methodology  is  ” worst”.  In  this  case  we  mean  worst  in  the  sense  that  difference  between  the 
Akaike  information  criterion  of  the  methodology  minus  the  null  method  is  maximal. 


22 


Figure  16:  Observed  spectrogram  (left)  and  predicted  spectrogram  (right)  of  the  five  patients  for  which 
the  proposed  methodology  does  ’’best” .  Where  we  mean  best  in  the  sense  that  the  AIC  of  the  proposed 
methodology  minus  AIC  of  the  null  method  is  minimal. 
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Figure  IT:  From  an  observed  spectrogram  (top)  the  methodology  returns  a  forward  model  designed  to 
best  the  data  (bottom  left).  Simply  omitting  the  contributions  to  the  predicted  spectrogram  from  aliasing 
and  tissue  motion  (bottom  right)  amounts  to  signal  processing  the  data  so  as  to  remove  these  artifacts. 


at  modeling  of  blood  flow,  because  models  generally  produce  flow  fields,  but 
TDU  observations  are  generally  not  of  the  flow  field  but  of  the  spectrogram. 
Expression  3,  therefore,  allows  one  to  verify  or  validate  models  of  blood  flow. 
For  instance,  closed-form  solutions  to  pulsatile  fluid  flow  in  an  elastic  tube 
date  from  the  1950’s  in  [2]  and  [1],  But,  more  recent  work  by  Myers[3] 
develops  analytic  solutions  for  a  curved  artery.  As  such,  one  could  test 
Myers’  model  by  simply  replacing  the  component  for  flow-field  calculation 
with  Myers’  solutions. 

Another  exercise  made  available  by  the  proposed  methodology  is  as  fol¬ 
lows:  All  of  the  data  sets  are  anonymized  but  it  appears  that  a  few  are  drawn 
from  repeated  patients;  that  they  are  drawn  from  a  single  subject  differing 
in  only  adjustment  of  the  TDU  device.  An  estimate  of  0  (t)  and  the  free 
parameters  of  the  spectrogram  calculation  are  made  as  usual  but  in  a  slight 
change  to  the  methodology  free  parameters  for  the  Windkessel  model  and  the 
flow  field  calculation  are  shared  between  the  two  data  sets.  In  other  words 
the  free  parameters  representing  the  physiology  of  the  patient  are  constant 
across  the  two  data  sets  and  optimized  for  the  data  sets  jointly  while  the 
free  parameters  representing  the  TDU  and  the  time  of  their  cardiac  cycles 
are  separate  and  optimized  for  each  data  set  independently.  The  resulting 
predicted  spectrograms  are  shown  in  the  middle  row  of  Figure  18,  and  the 
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corresponding  observed  spectrograms  shown  in  the  top  row.  Also,  as  long  as 
we  have  a  forward  model  representative  of  a  patients  physiology  we  may  do 
away  with  the  free  parameters  of  the  spectrogram  calculation  entirely  and 
replace  them  with  idealized  values.  In  this  application  the  flow  velocity  field 
could  be  thought  of  as  the  more  fundamental  quantity  and  the  simulated 
spectrogram  as  merely  a  flow  visualization  graphic.  This  is  shown  in  the 
bottom  row  of  Figure  18  Beyond  the  predicted  spectrograms  we  can  report 


Figure  18:  Two  data  sets  presumed  to  originate  from  the  same  patient  are  shown  at  the  top.  Shown  in 
the  middle  are  the  predicted  spectrogram’s  with  shared  Windkessel  models  and  flow  field  calculations.  In 
the  bottom  row  the  forward  models  are  identical  to  those  shown  in  the  middle  row,  with  the  exception 
that  the  spectrogram  calculation  free  parameters  values  are  replaced  with  idealized  values. 


(Table  4)  the  optimized  values  of  the  free  parameters  of  the  spectrogram 
calculation  corresponding  to  Figure  18.  With  the  caveat  that  we  do  not  ex¬ 
pect  to  have  found  a  global  optimum  for  these  free  parameters,  we  can  say 
that  out  of  all  the  sets  of  trial  parameter  values  these  are  the  most  likely. 
Without  a  methodology  like  the  one  of  this  paper  no  such  an  analysis  would 
be  possible. 
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parameter 


left 


right 


idealized 


focus  offset  x 
focus  offset  y 
focus  width 
off-bore 
aliasing  integer 
wall  spread 


0.003565  cm 
0.001167  cm 
0.003518  cm 
0.125520  rad 
2, 

6.904307  cm/s 


0.003163  cm 
0.002609  cm 
0.002170  cm 
0.342641  rad 
2, 

7.930820  cm/s 


0.0 

0.0 

a0 
0.0 
N/A 
7.21875  cm/s 


Table  4 


The  motivation  for  this  work  is  that  the  optimized  values  of  the  free  pa¬ 
rameters  will  be  found  useful,  but  optimization  and  attainment  of  their  most 
likely  values  remains  an  aim  for  future  work.  Fixed  parameter  values  which 
could  reliably  be  fit  to  data  would  have  novel  applications.  For  example,  al¬ 
though  not  explicitly  clear,  the  methodology  developed  here  leads  to  a  model 
of  autoregulation,  the  response  of  smooth  muscle  to  the  CO2  concentration 
in  blood  [25].  In  current  clinical  practice  the  visual  assesment  of  a  spec¬ 
trogram  is  an  important  tool  for  assessing  autoregulation.  In  the  proposed 
forward  model,  the  Young’s  modulus,  Poisson’s  ratio,  and  the  rest  diameter 
for  the  middle  cerebral  artery,  are  all  parameters  describing  the  properties 
and  state  of  the  smooth  muscle.  With  estimates  of  these  parameters,  the 
resulting  model  is  implicitly  a  model  of  autoregulation. 
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ORIGINAL  RESEARCH 


Mixture  Models  for  Estimating  Maximum 
Blood  Flow  Velocity 


Caren  Marzban,  PhD,  Wenxiao  Gu,  MS,  Pierre  D.  Mourad,  PhD 


Objectives — Agaussian  mixture  model  (GMM)  was  recently  developed  for  estimating 
the  probability  density  function  of  blood  flow  velocity  measured  with  transcranial 
Doppler  ultrasound  data.  In  turn,  the  quantiles  of  the  probability  density  function  allow 
one  to  construct  estimators  of  the  “maximum”  blood  flow  velocity.  However,  GMMs 
assume  gaussianity,  a  feature  that  is  not  omnipresent  in  observed  data.  The  objective  of 
this  work  was  to  develop  mixture  models  that  do  not  invoke  the  gaussian  assumption. 

Methods — Here,  GMMs  were  extended  to  a  skewed  GMM  and  a  nongaussian  kernel 
mixture  model.  All  models  were  developed  on  data  from  59  patients  with  closed  head 
injuries  from  multiple  hospitals  in  the  United  States,  with  ages  ranging  from  13  to  81 
years  and  Glasgow  Coma  Scale  scores  ranging  from  3  to  1 1.  The  models  were  assessed 
in  terms  of  the  log  likelihood  (a  goodness-of-fit  measure)  and  via  visual  comparison 
with  the  underlying  spectrograms. 

Results — Among  the  models  examined,  the  skewed  GMM  showed  a  significandy  (P  <  .05) 
higher  log  likelihood  for  56  of  the  59  patients  and  produced  maximum  flow  velocity 
estimates  consistent  with  the  observed  spectrograms  for  all  patients.  Kernel  mixture 
models  are  generally  less  “robust”  in  that  their  quality  is  inconsistent  across  patients. 
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Conclusions — Among  the  models  examined,  it  was  found  that  the  skewed  GMM  pro¬ 
vided  a  better  model  of  the  data  both  in  terms  of  the  quality  of  the  fit  and  in  terms  ofvisual 
comparison  of  the  underlying  spectrogram  and  the  estimated  maximum  blood  flow 
velocity.  Nongaussian  mixture  models  have  potential  for  even  higher-quality  assess¬ 
ment  of  blood  flow,  but  further  development  is  called  for. 

Key  Words — blood  flow;  brain;  head  injury;  noninvasive;  transcranial  Doppler  ultrasound 


The  importance  of  estimating  blood  flow  in  major  cerebral 
arteries  has  been  well  documented.1-5  The  reduction  of 
blood  flow  through  a  major  cerebral  artery,  such  as  caused 
by  cerebral  vasospasm,  can  lead  to  a  wide  range  of  disorders;  there¬ 
fore,  monitoring  blood  flow  has  important  clinical  consequences.6 
Transcranial  Doppler  ultrasound  imaging  is  one  of  the  methods  for 
assaying  blood  flow.7,8  Vasospasm  leads  to  reduced  blood  flow  but  an 
increased  blood  flow  velocity;  therefore,  estimating  the  maximum 
blood  flow  velocity  is  of  particular  importance.9 

The  time  series  of  a  flow  velocity  histogram  is  called  a  spectro¬ 
gram,  and  the  time  series  of  the  maximum  flow  velocity  is  called  an 
envelope.  Despite  their  clinical  importance,  neither  the  spectrogram 
nor  the  maximum  flow  velocity  can  be  observed  directly.  From  a  sta¬ 
tistical  perspective,  they  must  be  considered  population  parameters 
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that  are  to  be  estimated  from  data.  Moreover;  the  very 
notion  of  maximum  flow  velocity  is  ambiguous  because  the¬ 
oretically  flow  velocity  has  no  upper  bound.  Many  envelope 
estimation  algorithms  have  been  proposed10-15;  one  that  is 
frequently16  used  is  the  modified  geometric  method 
(MGM).10  That  method  essentially  searches  for  a  “kink”  in 
the  cumulative  histogram  of  flow  velocity  and  defines  the 
maximum  flow  velocity  as  the  flow  velocity  at  which 
the  kink  occurs.  It  provides  reliable  envelope  estimates  but 
lacks  the  important  feature  of  allowing  one  to  quantify  the 
uncertainty  in  maximum  flow  velocity.  By  contrast;  the  mix¬ 
ture  model  described  by  Marzban  et  al15  allows  one  to 
quantify  flow  velocity  in  terms  of  percentiles.  For  example; 
one  may  examine  the  envelope  corresponding  to  the  90th; 
95th;  and  99th  percentiles  of  flow  velocity.  Each  of  these 
provides  a  unique  and  well-defined  measure  of  “maximum” 
flow  velocity;  and  the  ability  to  do  this  is  the  main  advan¬ 
tage  of  mixture  models.  Furthermore;  the  maximum  flow 
velocity  as  defined  by  the  MGM  is  in  close  agreement  with 
that  produced  by  the  90th  percentile  envelope  in  a  mixture 
model.15  As  such;  mixture  models  enjoy  two  important 
features:  (l)  they  produce  envelopes  corresponding  to 
multiple  percentiles  of  flow  velocity;  and  (2)  one  of  the  per¬ 
centiles  (ie;  the  90th)  produces  an  envelope  that  is  approx¬ 
imately  equal  to  the  envelope  estimated  by  the  MGM. 

However;  the  specific  type  of  mixture  model  devel¬ 
oped  by  Marzban  et  al15  assumes  that  the  distribution  of 
flow  velocity;  at  any  point  in  time;  as  measured  by  a  tran- 
scranial  ultrasound  device;  is  gaussian.  The  agreement 
between  an  envelope  produced  by  such  a  gaussian  mixture 
model  (GMM)  and  that  produced  by  the  MGM  suggests 
that  either  the  distribution  of  flow  velocity  truly  is  gaussian; 
or  that  the  envelope  is  not  sensitive  to  violations  of  that 
assumption.  As  shown  below;  the  distribution  of  flow  veloc¬ 
ity  varies  with  time;  at  certain  times;  it  is  gaussian  (or  at  least; 
near  gaussian);  while  at  other  times  it  is  not.  The  deviations 
from  normality  manifest  themselves  mostly  as  a  skew. 
For  this  reason;  one  of  the  models  described  here  assumes 
that  the  underlying  distributions  of  flow  velocity  are  skewed 
gaussian.  Data  are  used  to  estimate  the  mean;  the  variance; 
and  the  shape  (or  skew)  parameter  of  the  distributions. 
Here;  this  model  is  referred  to  as  a  skewed  GMM.17-19 
Another  model  discussed  (briefly)  does  not  assume  a  para¬ 
metric  form  for  the  distribution  of  flow  velocity  at  all; 
instead;  the  distribution  is  inferred  by  using  kernel  esti¬ 
mates.20-24  Such  kernel  mixture  models  (KMMs)  are  flex¬ 
ible  in  that  they  accommodate  a  wide  range  of  distributions 
(including  gaussian).  However;  kernel  methods  involve  a 
smoothing  parameter  whose  value  is  a  priori  unknown. 
Generally;  extreme  values  of  the  parameter  correspond  to 


either  a  smooth  distribution  (eg;  gaussian)  or  a  highly 
irregular  distribution.  One  often  uses  some  criterion  (such 
as  the  maximization  of  likelihood)  to  optimize  the  smooth¬ 
ing  parameter. 

In  this  study;  the  skewed  GMM  and  KMM  were 
developed  for  estimating  the  envelope  for  each  of  59 
patients  in  our  data  set.  Although  quantitative  measures 
for  goodness  of  fit  were  used  to  compare  the  various  mod¬ 
els;  qualitative  comparisons  were  also  used.  This  practice  is 
necessary  because;  as  the  true  envelope  is  not  directly 
observable;  the  estimated  envelope  is  assessed  by  visually 
comparing  it  with  the  underlying  spectrogram. 

Materials  and  Methods 

Data 

The  data  for  this  work  were  collected  from  a  preclinical 
study  involving  59  patients  with  closed  head  injuries  from 
multiple  hospitals  in  the  United  States;  with  ages  ranging 
from  13  to  81  years  and  Glasgow  Coma  Scale  scores  rang¬ 
ing  from  3  to  1 1.  A  PowerMode  100  transcranial  Doppler 
device  (Spencer  Technologies;  Northborough;  MA)  was 
used  at  the  Harborview  Medical  Center  (Seattle;  WA); 
Columbia  University  Medical  College  (New  York;  NY);  and 
the  University  of  Texas  Medical  School  (Houston;  TX). 
The  device  has  an  ultrasound  carrier  frequency  that  ranges 
between  1.75  and  2  MHz  and  pulse  repetition  frequencies 
of 8000  to  1 0;000  Hz.  The  transcranial  Doppler  probe  was 
placed  on  the  temple  of  the  patient  by  the  sonographer 
after  application  of  ultrasound  gel.  In  all  cases;  the  sonogra¬ 
pher  sought  to  produce  an  optimal  spectrographic  signal; 
defined  here  as  one  with  a  maximum  systolic  flow  velocity 
in  the  main  branch  of  the  middle  cerebral  artery.  To  do 
SO;  the  sonographer  manipulated  the  angle  ofinsonation  of 
the  middle  cerebral  artery  by  the  probe;  the  power  of  the 
ultrasound;  as  well  as  the  depth  within  the  brain  from 
which  the  Doppler  spectral  data  were  obtained.  Further 
details  on  this  data  set  can  be  found  in  previous  studies  by 
Marzban  et  al.15’25  The  former  used  only  the  ultrasound- 
derived  transcranial  Doppler  spectra  for  the  purpose  of 
developing  the  aforementioned  GMMs  for  envelope  esti¬ 
mation.  In  the  latter;  the  data  were  used  to  develop  a  model 
of  arterial  blood  pressure;  which  can  in  turn  be  used  to  pre¬ 
dict  intracranial  pressure.  In  accordance  with  the  Institu¬ 
tional  Review  Board  for  each  hospital;  informed  consent 
was  obtained  from  all  patients  or  their  families. 

With  clinically  approved  transcranial  Doppler  units; 
blood  flow  velocity  in  the  middle  cerebral  artery  was  meas¬ 
ured  with  Doppler  ultrasound.  Data  were  collected  for 
periods  from  5  to  30  minutes.  All  retrospective  data  pro- 
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cessing  and  analysis  were  conducted  at  the  Applied  Physics 
Laboratory  of  the  University  of  Washington.  Although  the 
Doppler  spectral  time  series  were  initially  sampled  at  1 25  Hz, 
they  were  down-sampled  at  40  Hz.  This  resolution  was  suf¬ 
ficient  to  resolve  each  patient’s  systolic  rise,  diachrotic 
notch,  and  diastolic  minimum.  For  each  patient,  envelopes 
were  produced  over  a  fixed  duration  (118  sampled  points, 
or  about  3  minutes).  This  duration  is  sufficiently  long  to 
include  several  cardiac  cycles  while  allowing  details  of  the 
envelopes  to  be  visually  evident. 

Data  Processing  and  Statistical  Analysis 
In  previous  work,  it  was  shown  that  the  GMM  is  a  natural 
extension  of  the  MGM  in  that  the  former  provides  for 
multiple  estimators  of  maximum  flow  velocity  based  on  per¬ 
centiles  of  the  probability  density  function  of  the  flow  veloc¬ 
ity.  It  was  found  that  the  envelope  corresponding  to  the 
90th  percentile  of  the  probability  density  function  is 
comparable  with  the  envelope  estimated  by  the  MGM.15 
As  such,  the  GMM  is  superior  (to  the  MGM)  in  that  it 
allows  one  to  estimate  the  highest  blood  velocities,  and  effec¬ 
tively  their  uncertainty,  through  multiple  percentiles  of  the 
probability  density  function  of  flow  velocity.  The  MGM  is 
included  in  the  analysis  here  only  for  the  purpose  of  provid¬ 
ing  a  comparison  with  the  GMM.  In  fact,  since  the  GMM 
outperforms  the  MGM,  the  skewed  GMM  and  KMM 
were  both  compared  with  the  GMM.  The  KMM  leads  to 
envelopes  that  are  moderately  superior  to  those  based 
on  the  skewed  GMM  but  not  for  all  patients;  for  this  rea¬ 
son,  the  KMM  is  discussed  only  briefly  and  only  for  the 
purpose  of  discussing  future  work. 

As  explained  in  Appendix  1,  in  a  GMM,  the  distribution 
(more  accurately,  the  probability  density  function)  of  flow 


velocity  for  each  patient  and  at  any  time  is  assumed  to  be 
a  linear  combination  of  two  gaussians.  The  coefficients, 
called  mixing  coefficients,  as  well  as  the  parameters  of  the 
gaussians  are  estimated  from  data.  In  the  study  by  Marzban 
et  al,15  it  was  shown  that  the  two  gaussians  represent  the 
“signal”  and  the  “background”  portions  of  the  distribution, 
respectively.  Those,  then,  allow  one  to  use  the  upper  per¬ 
centiles  of  the  former  to  quantify  the  maximum  flow  velocity. 
The  left  panel  of  Figure  1  shows  an  instance  of  the  signal 
gaussian  (blue),  the  background  gaussian  (red),  the  90th, 
95th,  and  99th  percentiles  of  the  former  (blue  vertical 
lines),  and  the  maximum  flow  velocity  according  to  the 
MGM  (black  vertical  line).  The  data  (black  circles)  are 
from  a  single  patient  and  at  a  given  time. 

The  generalization  from  a  GMM  to  a  skewed  GMM 
is  relatively  straightforward  in  that  each  gaussian  in  the  mix¬ 
ture  model  is  replaced  by  a  skewed  gaussian  (Appendix  1 ) . 
This  process  introduces  only  1  more  parameter:  namely, 
the  shape  parameter,  which  affects  the  skew  of  the  gaussian. 
The  difference  in  the  models  may  be  visualized  by  the 
example  in  Figure  1.  The  middle  panel  of  Figure  1  shows 
the  various  components.  The  data  are  for  the  same  patient 
and  same  time  as  in  the  left  panel.  Note  the  evident  skew  in 
the  signal  and  the  background  components.  The  skewed 
GMM  algorithm  is  implemented  by  an  R  package  called 
mixsmsn,18'19  available  online  from  the  Comprehensive  R 
Archive  Network.26 

In  the  KMM  the  distribution  offlow  velocity  is  assumed 
to  be  a  linear  combination  of  two  distributions,  but  unlike 
the  GMM  and  skewed  GMM,  each  distribution  is  estimated 
by  using  kernel  methods  (Appendix  1 ).  The  kernel  method 
followed  here  uses  the  expectation  maximization  algorithm 
for  fitting  multivariate  nonpar ametric  mixtures  with  com- 


Figure  1.  Distribution  (more  accurately,  probability  density  function)  for  the  signal  (blue)  and  background  (red),  as  determined  by  the  GMM,  skewed 
GMM,  and  KMM.  The  blue  vertical  lines  denote  the  90th,  95th,  and  99th  percentiles  of  the  signal  flow  velocity  (FV),  and  the  black  vertical  line  marks 
the  maximum  flow  velocity  according  to  the  MGM. 
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pletely  unspecified  component  densities/9-21  except  for  a 
conditional  independence  assumption  described  by 
Benaglia  et  al.23  The  kernel  is  the  standard  normal  density 
function.  The  algorithm  is  implemented  in  an  R  package 
called  mixtools/3-25  available  from  the  Comprehensive  R 
Archive  Network. 

The  right  panel  of  Figure  1  shows  the  resulting  distri¬ 
butions.  Note  that  the  distributions  are  now  more  irregu¬ 
lar,  capturing  the  nonsmooth  structure  evident  in  the  data 
(black  circles).  Appendix  1  explains  the  manner  in  which 
the  smoothing  parameter  is  set. 

Results 

As  mentioned  previously,,  the  comparison  of  the  various 
models  examined  here  has  both  a  quantitative  component 
as  well  as  a  qualitative  (visual)  component.  The  "visual 
envelope”  suggested  by  a  visual  inspection  of  the  spectro¬ 
gram  is  contingent  on  the  color  coding  underlying  the 
spectrogram.15  Given  that  the  true  envelope  is  unknown, 
the  visual  envelope  can  be  considered  yet  another  estimate, 
and  despite  its  qualitative  nature,  the  visual  envelope  pro¬ 
vides  a  background  against  which  all  of  the  other  estimates 
may  be  viewed.  Then  a  model  is  declared  as  useful  if  the 
estimated  envelope  is  visually  consistent  with  the  underlying 
spectrogram.  Another  qualitative  criterion  that  is  applied  to 
assess  the  usefulness  of  the  models  is  the  consistency  (in 
time)  with  which  the  estimated  envelope  agrees  with  the 
spectrogram.  For  example,  if  the  estimated  envelope  has 
stand-alone  peaks  or  troughs  interspersed  throughout  the 
time  window  of  the  spectrogram,  then  the  model  is  deemed 
not  useful,  or  at  least  worthy  of  further  research. 

The  spectrograms  and  envelopes  for  the  GMM, 
skewed  GMM,  and  KMM  for  a  single  patient  are  shown 
in  Figure  2.  (The  corresponding  patient  is  different  than 


that  for  Figure  1;  the  reason  why  different  patients  were 
used  for  different  figures  was  to  illustrate  different  charac¬ 
teristics  of  the  various  methods.)  For  this  particular  patient, 
it  can  be  seen  that  the  3  percentiles  according  to  the  GMM 
(left  panel)  are  close  to  one  another  and  generally  lower  than 
the  envelope  based  on  a  visual  inspection  of  the  spectrogram. 
By  contrast,  the  3  envelopes  based  on  the  skewed  GMM 
are  more  widely  spread  and  more  consistent  with  the  spec¬ 
trogram  (middle  panel).  The  KMM’s  envelopes  (right 
panel)  have  similarities  with  and  differences  from  the  GMM 
and  the  skewed  GMM.  For  example,  they  underestimate 
the  flow  velocities  (as  for  the  GMM)  but  are  relatively 
spread  out  (as  for  the  skewed  GMM).  These  patterns  are 
generally  true  across  many  patients.  Also,  note  the  complete 
disagreement  between  the  spectrogram  and  the  99th  per¬ 
centile  envelope  at  the  time  just  before  20  (0.025  seconds). 
This  type  of  failure  is  characteristic  of  the  KMM,  and  it  is 
discussed  further  in  the  "Discussion”  section. 

Another  noteworthy  feature  in  Figure  2  is  that  GMM 
envelopes  are  generally  lower  than  the  aforementioned 
visual  envelope  but  only  for  relatively  high  flow  velocity 
values.  This  characteristic  is  evident  in  the  manner  in  which 
the  envelope  is  "below”  the  highest  flow  velocities  in  the 
spectrogram  but  only  when  flow  velocities  are  highest. 
The  skewed  GMM  solves  that  problem.  It  is  a  consequence 
of  the  fact  that  low  flow  velocity  values  are  consistent  with 
gaussian  distributions,  but  for  high  flow  velocity  values,  the 
distributions  are  more  skewed. 

The  time  dependence  of  the  shape  of  the  distributions 
can  be  confirmed  by  examining  the  values  of  the  shape 
parameter  as  a  function  of  time.  Figure  3  shows  the  time 
series  for  the  estimated  shape  parameter  for  the  signal  (left 
panel)  and  background  (right  panel).  It  can  be  seen  that 
for  most  of  the  time  series,  the  values  of  the  shape  param¬ 
eters  are  approximately  0;  ie,  the  distribution  of  the  flow 


Figure  2.  Spectrograms  (color  background)  and  90th,  95th,  and  99th  percentile  envelopes  (black)  according  to  the  GMM,  skewed  GMM,  and 
KMM  for  a  single  patient. 
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velocity  is  gaussian.  However,  there  are  also  times  at  which 
the  shape  parameters  are  unambiguously  not  0.  The  times 
correspond  to  regions  of  the  spectrogram  where  flow  veloc¬ 
ities  take  their  highest  values.  Note  that  when  the  shape 
parameters  are  not  0,  they  are  positive  for  the  signal  and 
negative  for  the  background  components  of  the  mixture 
model.  (A  positive  shape  parameter  value  corresponds  to 
a  right-skewed  gaussian.) 

As  seen  for  the  patient  corresponding  to  Figure  3,  the 
shape  parameter  is  0  for  most  times.  However,  that  is  not 
true  for  all  patients;  Figure  4  shows  a  histogram  of  the  esti¬ 
mated  shape  parameters  for  a  different  patient.  Note  the 
bimodal  structure;  the  left  mode  corresponds  to  the  back¬ 
ground,  whereas  the  right  mode  is  for  the  signal.  It  follows 
that  the  background  component  is  consistently  left  skewed 
(with  shape  parameters  in  the  -30  range),  and  the  signal 
component  is  either  gaussian  (corresponding  to  the  peak 
at  0)  or  right  skewed  (with  shape  parameters  between  0 
and  +30). 

Figures  such  as  those  in  Figure  2  are  useful  for  assess¬ 
ing  the  quality  of  an  envelope.  However,  they  were  imprac¬ 
tical  for  addressing  the  quality  of  the  envelopes  for  all  59 
patients.  Consequently,  although  we  performed  this  visual 
inspection  for  all  patients,  the  results  are  not  shown  here. 
Instead,  one  can  rely  on  the  goodness-of-fit  measure  itself, 
not  the  envelope,  to  compare  the  various  models  devel¬ 
oped  here.  There  are  numerous  quantities  that  can  be  used 
for  measuring  how  well  a  model  fits  the  data.  One  com¬ 
monly  used  measure  is  the  log  likelihood  (LL) .  Given  that 
the  fitting  of  the  model  is  performed  at  each  time,  there  exists 
a  time  series  for  the  LL  as  well  (not  shown) .  For  the  purpose 


of  comparing  two  models,  however,  it  is  sufficient  to 
examine  the  difference  between  LL  values  of  the  models. 
Therefore,  for  each  patient,  the  histogram,  or  box  plot,  of 
this  quantity  provides  a  visual  assessment  of  the  relative 
performance  of  the  models.  If  two  models  are  comparable 
in  terms  of  goodness  of  fit,  then  the  box  plot  will  be  centered 
on  the  number  0.  Otherwise,  the  models  can  be  deemed 
different.  It  is  important  to  point  out  that  the  LL  assesses 
the  goodness  of  fit  of  the  mixture  model  used  for  estimating 
an  envelope,  not  the  quality  of  the  estimated  envelope 
itself;  an  assessment  of  the  quality  of  the  estimated 
envelope  would  require  the  true  envelope,  which  is  unob¬ 
servable.  The  reason  for  evaluating  the  models  here  in 
terms  of  LL  and  the  quality  of  envelopes  is  to  allow  for  the 
possibility  that  one  model  may  outperform  another  model 
in  terms  of  LL  but  not  necessarily  in  terms  of  the  quality  of 
the  envelope;  or  vice  versa. 

Figure  5  shows  a  box  plot  of  LL  (skewed  GMM) 
minus  LL  (GMM)  for  all  59  patients.  (Recall  that  the  vari¬ 
ability  of  these  box  plots  is  due  to  time.)  The  fact  that  most 
box  plots  are  4  above”  the  vertical  line  at  0  implies  that  for 
most  of  the  patients,  the  LL  of  the  skewed  GMM  is  generally 
higher  than  that  of  the  GMM  across  time.  In  other  words, 
on  average  (across  time),  the  skewed  GMM  is  a  better 
model  of  the  data  than  the  GMM.  A  few  exceptions  are 
patients  19, 41,  and  53,  for  whom  the  box  plots  are  mostly 
centered  around  0,  and  as  such,  there  is  no  significant  dif¬ 
ference  between  the  GMM  and  skewed  GMM.  To  assess 
whether  the  skewed  GMM  has  higher  LL  values  than  the 
GMM,  a  1 -sided  t  test  was  performed.  All  of  the  P  values 
were  found  to  be  less  than  .05,  except  those  for  patients  19, 


Figure  3.  Time  series  of  the  estimated  shape  parameter  for  the  signal  and  background  distribution  for  a  single  patient. 

Signal  Background 
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41 ,  and  53  (shaded  box  plots  in  Figure  5).  (The  assump¬ 
tion  of  normality  for  the  t  test  is  confirmed  by  a  visual 
inspection  of  the  normal  quantile-quantile  plot.26) 

The  comparison  between  the  KMM  and  GMM  can 
be  done  in  the  same  fashion.  A  box  plot  for  the  difference 
between  their  LL  is  shown  in  Figure  6.  Given  the  systematic 
presence  of  the  boxes  to  the  right  of  the  vertical  line  at  LL  = 
0;  it  is  evident  that  the  KMM  has  a  generally  higher 
LL  than  the  GMM.  (All  of  the  P  values  were  <.05. )  It  is 
important  to  emphasize  that  this  conclusion  is  not  unex¬ 
pected.  The  KMM  allows  for  more  flexible  distributions 
and;  consequently;  leads  to  higher  LL  values  than  the 
GMM.  The  more  important  question  is  whether  the 
KMM  also  leads  to  better  envelopes  than  the  GMM;  and 
the  answer  to  that  is  no.  The  KMM  envelopes  for  all 
59  patients  were  compared  with  the  corresponding 
envelopes  from  the  GMM  (both  skewed  and  not)  in  terms 
of  a  visual  comparison  of  the  envelopes;  a  scatterplot  of  two 
envelopes  (eg;  Figure  5  of  Marzban  et  al15);  and  a  com¬ 
parison  of  scalar  summary  measures  (eg;  Figure  6  of 
Marzban  et  al15).  Although  the  details  of  the  comparison 
are  not  shown  here;  the  results  suggest  that  the  KMM 
envelopes  are  no  better  than  those  of  the  GMM  and  are 
often  comparable  with  those  produced  by  the  skewed 
GMM.  The  KMM  is  also  not  as  consistent  (across  time)  as 
the  GMM  or  skewed  GMM;  as  seen  by  the  “ peak”  in  the 
right  panel  of  Figure  2.  Therefore;  although  there  is  suffi¬ 
cient  evidence  to  advocate  the  use  of  the  skewed  GMM; 
that  is  not  the  case  for  the  KMM.  Further  discussion  of  the 
KMM  is  provided  in  the  next  section. 


Figure  4.  Histogram  of  the  estimated  shape  parameter  for  a  single 
patient.  The  left  (right)  mode  in  this  histogram  corresponds  to  the  back¬ 
ground  (signal)  component  of  the  mixture  model. 
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Discussion 

In  an  earlier  study;  it  was  shown  that  the  GMM  has  numer¬ 
ous  advantages  over  the  MGM  for  estimating  the  envelope. 
Here;  it  has  been  shown  that  a  skewed  GMM  has  the  same 
qualitative  advantages  as  the  GMM  (eg;  allowing  for  a 
percentile-based  notion  of  maximum  flow  velocity)  but 
also  outperforms  it  in  terms  of  both  the  quality  of  the  enve¬ 
lope  (assessed  through  a  visual  comparison  of  the  spec¬ 
trogram  and  the  envelopes)  and  a  quantitative  comparison 
(in  terms  of  the  goodness  of  fit).  One  may  wonder  whether 
the  higher  goodness-of-fit  values  for  the  skewed  GMM 
(compared  with  the  GMM)  may  be  a  consequence  of 
overfitting.  This  concern  is  unwarranted  because  the 
skewed  GMM  has  only  2  additional  parameters:  the  shape 
parameters  for  the  signal  and  background  components. 
As  such;  the  total  number  of  independent  parameters  to 
be  estimated  is  only  7  (mean;  standard  deviation;  and 
shape  parameter  for  each  of  the  signal  and  background 
components;  plus  1  for  the  mixing  coefficient).  There 
exists  a  large  body  of  theoretical  work  on  the  relationship 
between  the  number  of  parameters  and  the  number  of 
observations;  eg;  see  the  Vapnik-Chernovenkis  dimen¬ 
sion  in  the  work  by  Hastie  et  al27  (page  210).  On  the  prac¬ 
tical  side;  one  guideline  that  is  often  followed  is  from  the 
classic  bookby  Ryan28  (page  20);  who  recommends  at  least  10 
times  as  many  observations  as  the  number  of  parameters. 
Here;  the  data  used  for  estimating  the  7  parameters  number 
in  the  hundreds;  therefore;  overfitting  is  not  a  concern. 

An  issue  that  arises  in  the  mixture  model  approach  to 
envelope  estimation  is  the  identifiability  of  the  two  com¬ 
ponents  with  signal  and  background.  More  specifically;  a 
mixture  model  with  two  components  will  identify  two  com¬ 
ponents  in  the  distribution  offlow  velocity  values;  however; 
an  additional  criterion  must  be  introduced  to  decide  which 
of  the  components  should  be  identified  as  the  signal.  This 
factor  is  important  because  it  is  the  percentiles  of  the  signal 
component  that  lead  to  the  envelope.  Visually;  the  sig¬ 
nal  component  can  be  identified  as  the  one  on  the  left;  but 
quantifying  what  is  meant  by  "on  the  left”  is  nontrivial.  The 
approach  adopted  here  for  identifying  the  signal  and  back¬ 
ground  components  of  the  mixture  model  is  outlined  in 
Appendix  2.  That  proposal  is  somewhat  ad  hoc;  so  it  would 
be  desirable  to  develop  a  more  statistically  sound  method 
for  identifying  the  signal  component  of  mixture  models. 

As  described  above;  the  smoothing  parameter  of  the 
KMM  determines  the  shape  of  the  kernel  estimate  of 
the  probability  density  function.  Here,  an  algorithmic  cri¬ 
terion  is  used  to  set  the  values  of  that  parameter;  but  the 
criterion  optimizes  the  goodness  of  fit  (ie;  LL).  As  such; 
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the  criterion  does  not  automatically  lead  to  improved 
envelopes.  An  instance  of  this  is  shown  in  Figure  2,  where 
the  KMM  envelope  displays  an  abrupt  peak  at  approxi¬ 
mately  time  15  (0.025  sec),  which  is  clearly  inconsistent 
with  the  observed  spectrogram.  This  is  one  of  the  main  rea¬ 
sons  why  the  KMM  is  not  wholeheartedly  advocated  here. 


Another  reason  is  more  computational:  the  procedure  for 
estimating  kernels  is  rather  computer  intensive  and,  there¬ 
fore,  slow.  In  fact,  for  this  study,  the  KMM  was  developed 
on  data  for  which  the  lowest  5%  of  flow  velocity  values 

Figure  6.  Same  as  Figure  5  but  for  LL  for  the  KMM  minus  LL  for  the  GMM. 


Figure  5.  Box  plot  of  LL  for  the  skewed  GMM  minus  LL  for  the  GMM  for 
all  59  patients.  The  shaded  boxes  indicate  P>  .05. 
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were  trimmed  (deleted).  This  process  substantially 
increases  computational  speed  but  introduces  an  ad  hoc 
parameter  (ie,  the  5%).  In  short,  more  work  is  required 
before  one  can  benefit  from  the  flexibility  the  of  KMM  as 
a  useful  envelope  estimation  method. 

Appendix  1:  Mixture  Models 

In  a  mixture  model,  one  generally  assumes  that  the  data 
are  generated  from  an  underlying  probability  density  func¬ 
tion  [f(x)]  that  can  be  written  as  a  linear  combination  of 
more  well-known  probability  density  functions.  For  exam¬ 
ple,  for  the  case  of  two  components  (as  considered  in  this 
study),  the  probability  density  function  is  written  as 

(1)  f(x)  =  'Zhfk(x), 

where  A*  are  called  mixing  coefficients  and  must  sum  to  1. 
For  GMMs,/fc(x)  is  given  by 

+(v> 

where  p  and  O2  are  mean  and  variance  parameters,  and 
4>(.)  is  the  standard  gaussian  probability  density  function. 
Consequently  such  a  mixture  model  has  5  independent 
parameters:  the  mean  and  variance  parameters  for  each  of 
the  two  gaussian  probability  density  functions  and  one  for 
the  mixing  weights,  which  must  be  estimated  from  data. 

In  a  skewed  GMM,  the  component  probability  den¬ 
sity  functions  are  assumed  to  be  given  by 

(2)  fk(x)  =  2(j)  (^2_H.)  $  (a*), 

where  ®(x)  is  the  gaussian  cumulative  probability  distri¬ 
bution,  and  a  is  called  the  shape  (or  skew)  parameter.19 
Large  positive  values  of  the  shape  parameter  lead  to  a  prob¬ 
ability  density  function  with  a  more  pronounced  tail  for 
large  values  o(x  (ie,  a  right-skewed  distribution),  whereas 
large  negative  values  of  a  have  the  same  effect  but  for  large 
negative  values  of  x  (ie,  a  left-skewed  distribution).  In  this 
application,  the  skewed  GMM,  therefore,  has  2  more 
parameters  than  the  5  parameters  in  the  GMM.  Estimated 
values  of  the  shape  parameter  are  shown  in  Figure  3,  and 
the  skewed  probability  density  functions  are  shown  in  the 
middle  panel  of  Figure  1 .  The  presence  of  the  shape  param¬ 
eter  allows  one  to  fit  a  more  general  class  of  data  than  can  be 
fit  with  GMMs. 

Even  more  general  data  can  be  fit  with  kernel  meth¬ 
ods,  in  which  the  component  probability  density  functions 
in  the  mixture  model  are  themselves  assumed  to  be  a  sum 


of  the  form 


where  xt  are  the  observed  data;  n  is  the  sample  size;  and  the 
function  K(.),  which  is  called  the  kernel,  can  be  any  function 
that  integrates  to  1;  commonly,  and  in  the  present  applica¬ 
tion,  K(.)  is  taken  to  be  the  standard  gaussian  probability 
density  function.27  The  additional  parameter  h ,  called  the 
bandwidth,  effectively  determines  the  smoothness  of  the 
probability  density  function;  a  small  value  of  h  will  lead  to 
a  highly  "spiky,”  nonsmooth  probability  density  function, 
whereas  larger  values  of  h  will  lead  to  smoother,  more 
gaussian-looking,  probability  density  functions.  For  this 
study,  its  value  was  determined  as  follows:  The  default 
choice  of  the  bandwidth  h  minimizes  a  quantity  called 
asymptotic  mean  integrated  squared  error29  and  is  given  by 

(4)  h  =  0.9n~l^smin  j^SD, 

where  SD  and  IQR  are  the  standard  deviation  and  inter¬ 
quartile  range  of  all  n  observations,  respectively.  It  has  been 
shown  that  this  estimate  is  somewhat  smaller  than  the  true 
optimal  value,24  so  the  value  chosen  here  is  taken  to  be 
larger  than  the  default  value,  which  in  turn  leads  to 
smoother  fits.  As  a  result,  at  each  time,  and  for  each  patient, 
h-  10  is  used,  unless  equation  4  gives  a  larger  value,  in 
which  case  the  larger  value  is  selected.  The  fits  are  visually 
inspected  to  ensure  that  the  chosen  value  of  the  bandwidth 
does  not  lead  to  overfitting. 

Appendix  2:  Identifying  the  Signal  Distribution 

As  described  above,  the  mixture  model  identifies  two 
components  in  the  distribution  of  flow  velocity  values. 
However,  the  task  of  identifying  which  of  the  components 
is  the  signal  (and  which  is  the  background),  is  nontrivial. 
Here,  the  following  3  criteria  are  used  to  derive  an  algo¬ 
rithm  for  the  identification  task:  The  signal  component  is 
more  likely  the  one  with  (l)  a  lower  mean,  (2)  a  higher 
peak,  and  (3)  a  larger  standard  deviation.  The  means,  peaks, 
and  standard  deviations  of  the  two  components  are 
denoted  as  (pi,  p2),  (pn  P2);  and  (sb  s2),  respectively. 
The  ratios  are  defined  as  a  =  pi/ p2;  b  =  pi/p2;  and  c  =  si / s2. 
As  such,  the  signal  component  is  more  likely  to  be  the  com¬ 
ponent  labeled  1  when  a  is  small  and  b  and  c  are  large. 
Therefore,  ifx=(l/a)xf?xcis  larger  than  1,  then  the 
component  labeled  1  should  be  identified  as  the  signal. 
Similarly,  x  <  1  indicates  that  the  component  labeled  1 
should  be  identified  as  the  background. 
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A  Double-Gaussian,  Percentile-Based 
Method  for  Estimating  Maximum 
Blood  Flow  Velocity 


Caren  Marzban,  PhD ,  Paul  R.  Illian ,  BS ,  David  Morison,  BS ,  Pierre  D.  Mourad,  PhD 


Objectives — Transcranial  Doppler  sonography  allows  for  the  estimation  of  blood  flow 
velocity,  whose  maximum  value,  especially  at  systole,  is  often  of  clinical  interest.  Given 
that  observed  values  of  flow  velocity  are  subject  to  noise,  a  useful  notion  of  “maximum” 
requires  a  criterion  for  separating  the  signal  from  the  noise.  All  commonly  used  criteria 
produce  a  point  estimate  (ie,  a  single  value)  of  maximum  flow  velocity  at  any  time  and 
therefore  convey  no  information  on  the  distribution  or  uncertainty  of  flow  velocity. 
This  limitation  has  clinical  consequences  especially  for  patients  in  vasospasm,  whose 
largest  flow  velocities  can  be  difficult  to  measure.  Therefore,  a  method  for  estimating 
flow  velocity  and  its  uncertainty  is  desirable. 

Methods — A  gaussian  mixture  model  is  used  to  separate  the  noise  from  the  signal  dis¬ 
tribution.  The  time  series  of  a  given  percentile  of  the  latter,  then,  provides  a  flow  veloc¬ 
ity  envelope.  This  means  of  estimating  the  flow  velocity  envelope  naturally  allows  for 
displaying  several  percentiles  (eg,  95th  and  99th),  thereby  conveying  uncertainty  in  the 
highest  flow  velocity. 

Results — Such  envelopes  were  computed  for  59  patients  and  were  shown  to  provide 
reasonable  and  useful  estimates  of  the  largest  flow  velocities  compared  to  a  standard 
algorithm.  Moreover,  we  found  that  the  commonly  used  envelope  was  generally  con¬ 
sistent  with  the  90th  percentile  of  the  signal  distribution  derived  via  the  gaussian  mix¬ 
ture  model. 
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Conclusions — Separating  the  observed  distribution  of  flow  velocity  into  a  noise  com¬ 
ponent  and  a  signal  component,  using  a  double-gaussian  mixture  model,  allows  for  the 
percentiles  of  the  latter  to  provide  meaningful  measures  of  the  largest  flow  velocities 
and  their  uncertainty. 

Key  Words — blood  flow;  brain;  head  injury;  noninvasive;  transcranial  Doppler  sonography 

Transcranially  derived  Doppler  measurements  of  blood  flow 
in  major  cerebral  arteries  have  found  many  clinical  applica¬ 
tions.1  In  addition  to  assaying  stroke  risk  due  to  sickle  cell 
disease,  dysfunction  of  cerebral  autoregulation,  and  a  patent  fora¬ 
men  ovale,  among  other  etiologies,  some  of  the  earliest  applica¬ 
tions  targeted  monitoring  for  vasospasm  after  subarachnoid 
hemorrhage.2-5  Cerebral  vasospasm,  the  transient  reduction  of 
the  diameter  of  1  or  more  major  cerebral  arteries,  can  lead  to 
reduced  blood  flow  into  the  brain  and  hence  cerebral  ischemia  and 
an  increasing  chance  for  neurologic  damage.1  Monitoring  for  the 
onset  of  vasospasm  remains  an  important  application  of  transcra¬ 
nial  Doppler  sonography,  with  more  30,000  patients  per  year  requir¬ 
ing  daily  monitoring  for  1  week  or  more.6  Adding  to  this  primarily 
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civilian  population  are  military  patients  with  head  injuries 
after  exposure  to  explosions,  typically  roadside  bombs, 
with  half  of  these  patients  having  cerebral  vasospasm.7 

Transcranial  Doppler  sonography  measures  the  dis¬ 
tribution  of  blood  flow  speeds  within  a  blood  vessel  toward 
or  away  from  the  transducer,  with  negligible  flow  speeds 
adjacent  to  the  blood  vessel  wall  and  maximum  flow 
speeds  near  the  center  of  the  vessel.8  Critical  for  the  assay 
of  cerebral  vasospasm,  among  other  uses  of  transcranial 
Doppler  sonography,  is  successful  capture  of  the  speed  of 
the  fastest  flowing  blood  within  the  major  cerebral  arter¬ 
ies,  since  this  value  acts  as  a  useful  proxy  for  blood  vessel 
narrowing.1  Capturing  this  speed  is  a  particularly  chal¬ 
lenging  problem,  since  vasospasm  reduces  the  volume  of 
blood  flow  while  accelerating  the  blood  flow  speed,9  hence 
reducing  the  target  for  transcranial  Doppler  measurements 
while  straining  against  the  upper  limits  of  ultrasound  data 
processing  due  to  aliasing.10 

The  time  series  of  a  flow  velocity  histogram  is  gener¬ 
ally  referred  to  as  a  spectrogram,  and  the  time  series  of  the 
maximum  flow  velocity  is  called  an  envelope.  Although  the 
spectrogram  conveys  a  great  deal  of  information,  the  enve¬ 
lope  is  often  the  only  quantity  a  clinician  examines.  This 
practice  is  reasonable  because  the  information  contained  in 
a  spectrogram  can  be  displayed  in  different  ways,  leading  to 
different  conclusions.  For  example,  Figure  1  shows  a  spec¬ 
trogram  with  2  different  color  schemes;  in  the  top  panel, 
the  spectrogram  is  linearly  related  to  the  color  scale, 
whereas  in  the  bottom  panel,  the  colors  are  proportional  to 
the  square  root  of  the  spectrogram. 

Figure  1.  Spectrograms  of  2  different  mappings  for  assigning  color  to 
the  flow  velocity  (FV)  values:  linear  (top)  and  sguare  root  (bottom). 


This  inherent  ambiguity  in  the  information  gleaned 
from  a  spectrogram  also  reflects  itself  in  the  corresponding 
envelope.  In  practice,  observed  flow  velocity  values  can 
range  from  1  to  300  cm/ s.  Therefore,  to  obtain  a  useful 
estimate  of  maximum  flow  velocity,  one  must  introduce 
some  criterion  that  defines  what  is  meant  by  maximum. 
Many  such  criteria  (standards  for  transcranial  Doppler 
measurements11-15)  are  based  on  the  cumulative  his¬ 
togram  of  flow  velocity.  Figure  2  shows  an  example  of  the 
relative  frequency  histogram  of  flow  velocity  (top)  and 
the  corresponding  cumulative  histogram  (bottom)  at  a 
specific  time  for  a  specific  patient  (hereafter,  patient  X); 
the  latter  is  obtained  by  integrating  the  former  from  left  to 
right.  One  may  define  the  maximum  flow  velocity  as  the 
value  at  which  the  cumulative  histogram  "levels  off,”  but 
there  exist  different  criteria  corresponding  to  different 
objective  measures  of  that  point.11-15  The  vertical  bar  in 
Figure  2  marks  the  flow  velocity  at  which  it  is  maximum 
according  to  the  modified  geometric  method.11'12  The  top 
panel  in  Figure  3  shows  the  spectrogram  and  the  modified 
geometric  method  envelope  (in  black)  for  patient  X. 

Given  that  the  proposed  algorithm  is  compared  with 
the  modified  geometric  method  algorithm,  a  brief  review 
of  the  latter  is  in  order.  The  modified  geometric  method 
algorithm  and  its  variants11'12  effectively  rotate  (clockwise 
and  about  the  origin)  the  cumulative  histogram  by  some 
amount.  The  effect  of  such  a  rotation  is  that  the  point  at 
which  the  cumulative  histogram  levels  off  translates  to  a 
point  at  which  the  rotated  cumulative  histogram  reaches 

Figure  2.  Histograms  of  flow  velocity  (FV)  for  patient  X  at  a  given  instant 
in  time  (top),  and  the  corresponding  cumulative  (ie,  integrated)  his¬ 
togram  (bottom). 
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the  maximum.  Given  certain  smoothness  constraints,  it  has 
been  shown  that  this  point  corresponds  to  the  maximum 
flow  velocity.12  Note  that  the  modified  geometric  method 
algorithm  generates  only  point  estimates  of  the  maximum 
flow  velocity  and  so  provides  no  natural  measure  of  uncer¬ 
tainty.  This  limitation  is  evident  in  the  cumulative  his¬ 
togram  in  Figure  2:  Although  the  modified  geometric 
method  provides  an  estimate  of  the  point  at  which  the 
cumulative  histogram  levels  off,  quantifying  the  corre¬ 
sponding  uncertainty  is  by  no  means  unambiguous.  Con¬ 
veying  uncertainty  is  important  because  different  levels  of 
uncertainty  call  for  different  clinical  actions. 

The  main  proposals  in  this  study  were  ( 1 )  to  develop 
a  new  means  of  separating  the  noise  and  the  signal  com¬ 
ponents  of  a  spectrogram  and  (2)  to  quantify  the  largest 
values  in  the  signal  using  the  concept  of  percentile.  The  nth 
percentile  of  a  quantity  x  is  defined  as  the  value  of  x  for 
which  n%  of  the  x  values  are  smaller.  For  example,  the  95  th 
percentile  offlow  velocity  is  the  flow  velocity  value  at  which 
95%  of  flow  velocity  values  are  smaller.  An  envelope  based 
on  the  95th  percentile  of  the  histogram,  then,  provides  a 
meaningful  estimate  of  the  “top  5%”  of  flow  velocities. 
Moreover,  2  (or  more)  envelopes  at  different  percentiles 
can  convey  some  notion  of  uncertainty.  In  this  study, 
envelopes  based  on  3  percentiles  were  considered:  90th, 
95th,  and  99th. 

Figure  3.  Top,  Spectrogram  for  patient  X  (color  background)  and  4  esti¬ 
mates  of  the  envelope:  modified  geometric  method  envelope  (black) 
and  90th,  95th,  and  99th  percentile-based  envelopes.  Bottom,  Modi¬ 
fied  geometric  method  envelope  (black),  95th  percentile-based  enve¬ 
lope  (red),  and  90th  and  99th  percentile-based  envelopes  (blue).  FV 
indicates  flow  velocity. 
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The  use  of  a  percentile  to  quantify  the  largest  flow 
velocity  presumes  that  the  histogram  is  the  correct  his¬ 
togram  offlow  velocity.  In  practice,  observed  spectrograms 
are  contaminated  by  various  types  of  noise.  One  of  these 
sources  of  noise  is  evident  in  the  small  “hump”  appearing 
on  the  right  side  of  the  (relative  frequency)  histogram  in 
Figure  2.  Flow  velocity  values  and  their  associated  per¬ 
centiles  are  useful  only  if  they  pertain  to  the  non-noise 
component  of  the  histogram,  hereafter  called  the  signal  dis¬ 
tribution.  In  fact,  the  90th  percentile  of  the  histogram 
shown  in  Figure  2  is  near  the  right-most  hump  and  there¬ 
fore  unrealistic.  Therefore,  to  compute  useful  percentile- 
based  envelopes,  one  must  first  disambiguate  the  signal 
and  noise  contributions  to  the  histogram.  To  that  end, 
here,  a  gaussian  mixture  model  with  2  components16,17  was 
used  to  represent  the  noise  and  signal  components,  respec¬ 
tively.  The  component  appearing  to  the  left  (closest  to  the 
origin)  was  defined  as  the  signal  distribution.  Armed  with 
the  signal  distribution,  the  3  percentiles  were  computed  at 
each  time,  and  their  time  series  was  computed  for  59 
patients.  The  3  percentile-based  envelopes  for  patient  X 
are  shown  in  Figure  3,  both  with  the  spectrogram  (top 
panel),  and  without  (bottom  panel). 

In  this  article,  the  details  of  this  percentile-based 
approach  for  estimating  an  envelope  are  presented.  It  was 
found  that  the  estimates  are  visually  consistent  with  the 
underlying  spectrograms,  and  the  modified  geometric 
method  envelope  is  approximately  consistent  with  the 
95th  percentile  envelope.  Also,  it  is  shown  that  percentile- 
based  envelopes  naturally  allow  for  displaying  envelope 
uncertainty. 

Materials  and  Methods 

Data 

Patient  data  for  this  preclinical  study  were  collected  from 
a  variety  of  hospitals  in  the  United  States,  following  a  study 
at  the  University  of  Washington  led  by  Dr  Mourad.  The 
data  included  arterial  blood  pressure  and  intracranial  pres¬ 
sure,  but  those  data  were  not  used  for  this  study;  only  sono- 
graphically  derived  transcranial  Doppler  spectra  were  used. 
Further  details  on  the  patient  data  maybe  found  in  an  arti¬ 
cle  by  Marzban  et  al.18  In  accordance  with  the  Institutional 
Review  Board  for  each  hospital,  informed  consent  was 
obtained  from  all  patients  or  their  families. 

Doppler  spectral  time  series  ofblood  flow  speed  in  the 
middle  cerebral  artery  (or  flow  velocity  -  flow  velocity) 
were  acquired  at  each  institution  via  clinically  approved 
transcranial  Doppler  units.  Data  acquisition  lengths  varied 
from  5  to  30  minutes.  All  retrospective  data  processing  and 
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analysis  were  conducted  at  the  Applied  Physics  Laboratory 
of  the  University  of  Washington.  Doppler  spectral  time 
series  data  collected  via  the  Applied  Physics  Laboratory's 
hospital  cohort  were  initially  digitized  at  125  Hz  in  time  and 
with  a  frequency  resolution  of  128  Hz  at  a  given  time. 
We  then  down-sampled  these  Doppler  time  series  to  a 
resolution  of  40  Hz  in  time,  sufficient  to  resolve  all  of  the 
following  features  in  all  of  our  patients'  data:  systolic  rise, 
diachrotic  notch,  and  diastolic  minimum.  We  also  selected 
a  fixed  duration  (118  time  steps,  or  ~3  minutes)  from  each 
of  the  59  patients  for  statistical  analysis.  This  duration  was 
sufficiently  long  to  include  several  cardiac  cycles  while  still 
allowing  details  of  the  envelopes  to  be  visually  evident. 


Data  Processing  and  Statistical  Analysis 
As  mentioned  in  the  introduction,  the  main  aims  of  the 
proposed  method  were  to  first  represent  the  instantaneous 
Doppler  frequency  distribution  of  the  flow  velocity  at  a 
given  moment  in  time  and  then  to  produce  spectral 
envelopes  of  those  frequency  distributions  through  time 
based  on  percentiles  of  that  distribution  in  the  Doppler  fre¬ 
quency.  Distinguishing  the  signal  from  the  noise  requires 
a  model.  In  this  article,  it  is  assumed  that  the  underlying 
distribution  offlow  velocity  consists  of  2  distributions,  cor¬ 
responding  to  the  signal  and  noise,  and  that  both  are  gauss- 
ian.  This  type  of  model  is  a  special  case  of  gaussian  mixture 
models,  wherein  a  distribution  is  assumed  to  be  a  linear 
combination  of  gaussian  distributions.16,17  The  weights 
in  the  linear  combination  (called  mixing  proportions)  and 
the  parameters  of  the  gaussian  distributions  are  then  esti¬ 
mated  from  data  via  some  optimization  procedure,  (in  this 
study,  the  expectation-maximization  algorithm  was  used, 
but  other  parameter  estimation  methods  are  equally  ade¬ 
quate.  The  expectation-maximization  algorithm  maxi¬ 
mizes  the  conditional  expected  log  likelihood.19)  Figure  4 
shows  the  signal  component  (black)  and  the  noise  com¬ 
ponent  (red)  for  the  example  shown  in  Figure  2.  Also 
shown  are  the  locations  of  the  maximum  flow  velocity 
according  to  modified  geometric  method  (black  vertical 
line)  and  the  90th,  95th,  and  99th  percentiles  of  the  signal 
distribution  (blue  vertical  lines).  All  statistical  analyses 
were  performed  in  R.20 

Results 


The  previous  section  illustrates  the  proposed  method  on  a 
given  patient  and  at  a  given  time.  The  top  panel  in  Figure 
3  shows  the  3  percentile-based  envelopes  (blue  curves)  for 
patient  X.  It  is  evident  that  the  percentile-based  envelopes 
are  consistent  with  the  underlying  spectrogram.  The  bot¬ 


tom  panel  in  Figure  3  shows  4  envelopes;  the  modified 
geometric  method  envelope  (black)  is  comparable  with  the 
95th  percentile  envelope  (red),  sandwiched  between 
the  90th  and  99th  percentile  envelopes  (blue).  Clearly, 
displaying  multiple  percentile-based  envelopes  conveys  a 
sense  of  the  uncertainty  in  the  maximum  flow  velocities. 

Although  technical,  it  is  also  worth  mentioning  2 
additional  steps  taken  to  generate  the  percentile-based 
envelopes  shown  in  Figure  3:  (l)  The  optimization  algo¬ 
rithm  used  for  the  gaussian  mixture  model  requires  initial 
estimates  of  the  parameters  of  interest.  At  time  0,  the  initial 
values  of  the  parameters  are  assigned  randomly,  and  the 
optimization  algorithm  produces  parameter  values  that 
characterize  the  2  best-fit  gaussian  components  of  the  flow 
velocity  histogram  at  that  time.  At  all  "future”  times,  the 
initial  values  of  the  parameters  are  set  to  the  final  values 
obtained  at  the  previous  time  step.  This  process  allows 
some  of  the  serial  (auto)  correlation  in  the  spectrogram  to 
be  reflected  in  the  envelopes.  If  the  values  of  the  gaussian 
mixture  model  parameters  are  initialized  randomly  for  every 
and  all  time  steps,  the  resulting  envelopes  are  somewhat 
less  smooth  than  those  shown  in  Figure  3.  (2)  A  running- 
median  filter  with  a  window  size  of  3  seconds  is  applied  to 
the  envelopes  to  smooth  them  even  further.  The  size  of  the 
window  (ie,  3)  is  obtained  by  trial  and  error  and  by  a  visual 
comparison  with  the  modified  geometric  method  enve¬ 
lope.  The  modified  geometric  method  envelope  too  is 
smoothed  by  a  variety  of  techniques,  eg,  an  averaging  (over 
multiple  cardiac  cycles)  filter  and  a  finite  impulse  response 


Figure  4.  Histogram  of  flow  velocity  (FV)  for  patient  X  (circles)  and  the 
estimated  densities  for  the  2  gaussian  components.  The  black  vertical 
bar  denotes  the  maximum  flow  velocity  according  to  the  modified  geo¬ 
metric  method,  and  the  blue  vertical  bars  correspond  to  the  90th,  95th, 
and  99th  percentiles  of  the  signal  distribution  on  the  left. 
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filter.  Although  the  smoothness  of  the  displayed  envelopes 
is  important  in  a  clinician  s  decision  making,  it  is  a  feature 
that  is  easily  controlled  (eg,  by  a  single  parameter,  such  as 
the  window  size  of  the  running-median  filter)  and  so  is  not 
of  serious  concern. 

The  procedure  was  applied  to  all  59  patients  in  the 
data  set.  It  is  impractical  to  show  all  59  figures  analogous  to 
Figure  3.  Instead,  the  3  percentile-based  envelopes  were 
compared  with  the  modified  geometric  method  envelope 
using  scatterplots,  which  were  then  further  summarized 
into  scalar  measures. 

Each  panel  in  Figure  5  shows  the  scatterplot  of  one  of 
the  percentile-based  envelopes  versus  the  modified  geo¬ 
metric  method  envelope  for  patient  X.  There  are  1 18  points 
on  each  scatterplot,  corresponding  to  the  118  time  steps 
displayed  in  Figure  3.  The  diagonal  line  is  simply  a  line  of 
slope  1,  y-intercept  0.  It  is  clear  from  the  middle  panel  that 
the  95th  percentile  envelope  generally  agrees  with  the  mod¬ 
ified  geometric  method  envelope.  As  expected,  the  90th  and 
99th  percentile  envelopes  are  generally  lower  and  higher, 
respectively,  than  the  modified  geometric  method  enve¬ 
lope.  The  disagreements  between  the  percentile-based 
envelopes  and  the  modified  geometric  method  envelope 
are  largest  for  higher  flow  velocity  values.  This  result  is  espe¬ 
cially  evident  in  the  99th  percentile  envelope  (bottom 
panel)  in  the  manner  in  which  the  deviation  of  the  scatter¬ 
plot  from  the  diagonal  line  is  largest  for  higher  flow  veloc¬ 
ity  values.  In  other  words,  the  systolic  peaks  of  the  99th 
percentile  envelope  are  in  fact  higher  than  might  be 
expected  from  the  modified  geometric  method  envelope. 

One  may  further  summarize  each  of  these  panels  by 
computing  the  root  mean  squared  error  between  the  mod¬ 
ified  geometric  method  envelope  and  each  percentile- 
based  envelope.  For  patient  X,  the  root  mean  squared  error 
values  corresponding  to  the  3  panels  in  Figure  5  are  17.9, 
16.9,  and  31.2  cm/s.  In  addition  to  confirming  that  the 
modified  geometric  method  envelope  is  closest  to  the  95th 
percentile  envelope  (for  patient  X),  these  numbers  have 
useful  interpretations  as  well;  eg,  it  can  be  said  that  the 
modified  geometric  method  envelope  and  the  95th  per¬ 
centile  envelope  have  a  typical  deviation  of  about  16.9 
cm/ s  across  the  118  time  steps  examined. 

Although  useful  for  comparing  different  envelopes, 
root  mean  squared  error  values  do  not  assess  whether  the 
difference  in  the  envelopes  is  due  to  the  amount  of  scatter 
in  the  scatterplot  or  to  an  overall  shift.  There  exists  a 
decomposition  of  root  mean  squared  error  that  allows 
one  to  examine  both  contributions.  The  details  of  the 
decomposition  are  presented  in  “Appendix:.”  Here,  suf¬ 
fice  it  to  say  that  the  correlation  coefficient  (denoted 


Corr)  and  bias  (ie,  mean  [modified  geometric  method]  - 
mean  [percentile-based  method] )  gauge  the  amount  of 
scatter  and  the  overall  shift,  respectively.  The  Corr  values 
for  the  3  panels  in  Figure  5  are  0.80,  0.81,  and  0.81  (from 
top  to  bottom),  suggesting  that  the  amount  of  scatter  in 
the  3  panels  is  comparable.  Said  differently,  any  of  the  3 
percentile-based  envelopes  provides  an  adequate  approx¬ 
imation  to  the  modified  geometric  method  envelope, 
if/when  correlation  is  the  measure  of  agreement.  The  main 
difference  between  the  3  panels  is  in  the  Bias  values:  8.5, 
-2.6,  and  -23.5  (from  top  to  bottom).  These  numbers  have 
useful  interpretations,  too;  eg,  one  can  say  that  the  95th 
percentile  envelope  is  generally  shifted  above  the  modified 
geometric  method  envelope  by  about  2.6  cm/ s  for  patient  X. 

To  compare  the  various  envelopes  across  all  of  the 
patients,  root  mean  squared  error,  Corr,  and  Bias  were 
computed  for  all  59  patients.  Their  distributions  are  shown 
in  Figure  6.  The  left-most  box  plot  in  the  top  panel  sum¬ 
marizes  the  histogram  (across  59  patients)  of  the  root 
mean  squared  error  computed  for  the  90th  percentile 
envelope  and  the  modified  geometric  method  envelope 
(across  118  time  steps). 

The  middle  panel  in  Figure  6  is  analogous  to  the  top 
panel  except  that  the  measure  of  agreement  between  the 
modified  geometric  method  envelope  and  the  percentile- 
based  envelope  is  Corr.  The  3  box  plots  are  comparable, 
implying  that  all  3  percentile-based  envelopes  are  compa¬ 
rably  correlated  with  the  modified  geometric  method 
envelope  across  all  59  patients.  The  bottom  panel  shows 
that,  on  average  (across  all  59  patients),  the  90th  and  99th 

Figure  5.  Scatterplots  of  the  3  percentile-based  envelopes  and  the  mod¬ 
ified  geometric  method  envelope  for  patient  X.  FV  indicates  flow  velocity. 
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percentile  envelopes  are  below  and  above  the  modified 
geometric  method  envelope,  respectively.  By  comparing 
the  3  box  plots  in  the  top  or  bottom  panel,  it  can  be  seen 
that  the  agreement  between  the  modified  geometric 
method  envelope  and  the  95th  percentile  envelope,  noted 
for  patient  X,  extends  to  all  59  patients. 

It  is  worth  displaying  the  envelope  for  a  patient  for 
whom  such  scalar  measures  are  generally  poor.  Figure  7 
shows  the  spectrogram  and  envelopes  for  the  patient  with 
the  lowest  Corr  value.  It  is  clear  that  the  modified  geo¬ 
metric  method  envelope  (black)  is  in  fact  incorrect  and 
that  the  percentile-based  envelopes  (blue)  are  more  con¬ 
sistent  with  the  underlying  spectrogram.  In  other  words, 
the  low  Corr  value  for  this  patient  does  not  imply  a  defect 
on  the  part  of  the  percentile-based  envelope  but  rather  a 
defect  in  the  modified  geometric  method  envelope. 

Discussion 

It  has  been  shown  that  a  double-gaussian  mixture  model 
can  be  used  to  identify  a  component  of  the  histogram 
of  flow  velocity  corresponding  to  a  "signal.”  In  turn, 
percentiles  of  this  signal  distribution  provide  envelopes, 
which  meaningfully  quantify  the  largest  flow  velocity. 
Not  only  are  these  percentile-based  envelopes  visually 
consistent  with  their  underlying  spectrograms,  they 
are  also  objectively  consistent  with  a  commonly  used 
envelope  based  on  the  modified  geometric  method.  The 
percentile-based  envelopes  not  only  objectively  quantify 


what  is  meant  by  "largest  flow  velocities”  but  also  have  the 
added  advantage  (over  existing  envelopes)  of  allowing 
for  a  natural  display  of  uncertainty  in  the  envelopes. 

Given  that  the  analysis  is  performed  on  real  (not  sim¬ 
ulated)  data,  the  true  envelope  is  not  known.  As  such,  the 
quality  of  each  envelope  cannot  be  assessed  objectively  on 
its  own  merit.  For  that  reason,  the  analysis  here  has  been 
based  on  visually  comparing  an  envelope  with  the  under¬ 
lying  spectrogram  (eg,  as  in  Figure  3),  or  objectively  com¬ 
paring  each  of  the  percentile-based  envelopes  with  the 
modified  geometric  method  envelope  (eg,  Figures  5  and 
6).  In  spite  of  the  central  role  played  by  the  modified 
geometric  method  envelope  in  comparing  envelopes,  it  is 
important  to  recall  that  the  modified  geometric  method 
envelope  is  not  the  true  envelope  (eg,  see  Figure  7).  Here, 
it  is  used  as  a  standard  only  because  of  its  common  usage. 

Of  the  3  percentile-based  envelopes,  the  95th  per¬ 
centile  envelope  is  closest  to  the  modified  geometric 
method  envelope,  but  the  agreement  is  not  perfect,  as  seen 
by  the  box  plots  in  Figure  6.  For  example,  according  to  the 
middle  box  plot  in  the  bottom  panel,  on  average  (across 
patients),  the  95th  percentile  envelope  is  slightly  above  the 
modified  geometric  method  envelope.  One  may  ask  what 
is  the  exact  percentile  corresponding  to  the  modified  geo¬ 
metric  method  envelope,  but  the  question  itself  is  inap¬ 
propriate  because  it  elevates  the  status  of  the  modified 
geometric  method  envelope  to  a  "reference  standard,” 
when  in  fact  it  is  not.  What  is  more  important  is  that  the 
proposed  approach  produces  a  useful  spectrogram  as  well 


Figure  6.  Distribution  of  the  root  mean  squared  error  (RMSE),  correla¬ 
tion  coefficient  (Corr),  and  Bias  (defined  in  “Appendix”)  for  all  59  patients 
examined  in  this  work. 


Figure  7.  Same  as  Figure  3,  but  for  a  patient  for  whom  the  modified  geo¬ 
metric  method  envelope  (black)  is  clearly  incorrect.  FV  indicates  flow 
velocity. 
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as  an  interpretable  (in  terms  of  percentiles)  generalization 
of  approaches  that  naturally  produce  only  a  single  enve¬ 
lope  (eg;  modified  geometric  method). 

A  percentile-based  envelope  is  based  on  a  percentile  of 
the  signal  distribution;  in  which  the  signal  (and  noise)  dis¬ 
tribution  is  defined  from  fitting  a  gaussian  mixture  model 
with  2  components  to  the  whole  histogram;  at  a  given  time. 
There  is  some  justification  for  examining  the  3  compo¬ 
nents.  For  example;  in  Figure  4,  one  may  argue  that  there 
are  possibly  3  overlaying  histograms:  1  corresponding  to 
the  large  hump  on  the  left  (ie;  the  signal);  1  associated  with 
the  small  hump  on  the  right  (ie;  associated  with  an  alias¬ 
ing  reflection  on  Doppler  imaging;  and  another  corre¬ 
sponding  to  a  uniform  "background”  spanning  the  full 
range  of  flow  velocity  values.  We  repeated  the  entire  analy¬ 
sis  but  with  3  gaussian  components.  The  results  are  incon¬ 
clusive  and  require  further  research.  For  some  patients;  the 
results  do  not  change  substantially;  but  for  others;  they  do 
in  ways  that  can  be  considered  "better”  or  "worse”  depend¬ 
ing  on  the  measure  of  quality. 

The  percentiles  of  the  signal  distribution  allow  for 
meaningful  envelopes;  which  can  objectively  assess  the 
largest  flow  velocity.  In  addition;  displaying  multiple 
envelopes  can  convey  information  on  the  uncertainty  in  the 
observed  flow  velocity.  One  can  envisage  an  alternative 
measure  of  uncertainty.  For  example;  one  may  consider  the 
envelope  corresponding  to  a  fixed  percentile;  say  95th;  and 
then  obtain  the  distribution  of  that  quantity  via  some 
resampling  method.21  In  turn;  that  (sampling)  distribution 
can  be  used  to  compute  confidence  intervals  for  the  true 
95th  percentile  envelope.  Such  interval  estimates  of  the 
envelope  can  be  useful  in  conveying  uncertainty  if/ when  a 
specific  percentile  is  of  interest.  Otherwise;  displaying  mul¬ 
tiple  envelopes  corresponding  to  different  percentiles;  as 
done  here;  is  sufficient  for  conveying  uncertainty. 

Appendix 

Let  x.  and  y.  denote  the  values  of  the  modified  geometric 
method  envelope  and  a  percentile-based  envelope;  respec¬ 
tively;  at  time  i.  The  root  mean  square  error  between  2 
envelopes  is  defined  as  the  square  root  of  the  mean  square 
error  ( MSE ) : 

(!) 

where  n  is  the  length  of  the  time  steps  (here;  1 18).  It  can  be 
shown  that 

(2)  MSE  =  s2  +  s2  -  2(Corr \sy  +  /L  (Bias)2, 


where  s  and  s  are  the  sample  standard  deviations  of  x  and 

x  y  -t 

y,  respectively;  Corr  denotes  the  Pearson  correlation  coef¬ 
ficient;  and 

(3)  Bias  =  mean  x  -  meany. 

According  to  Equation  a  comparison  of  2  envelopes  in 
terms  of  mean  square  error  (or  root  mean  square  error)  is 
equivalent  to  a  comparison  in  terms  of  Corr  and  in  terms 
of  Bias;  the  former  compares  the  correlation  between  the 
envelopes;  and  the  latter  measures  the  overall  shift  between 
them. 
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